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Page  114       Symbols  used  in  figure  Fl  are: 

PEXIT    =  refrigerant  pressure  at  tube  exit 
PFLASH  =  refrigerant  flash  pressure 

PO         =  refrigerant  pressure  before  expansion  device  inlet 
PI         =  refrigerant  pressure  in  expansion  device  inlet 
X  =  refrigerant  quality  before  expansion  device  inlet 

Page  116       Symbols  used  in  figure  Gl  are: 

H      =  refrigerant  enthalpy  at  tube  end  obtained  from  present 
calculation  loop 
=  refrigerant  enthalpy  at  tube  end  obtained  from  previous 

calculation  loop 
=  refrigerant  enthalpy  at  coil  exit  obtained  from  present 

calculation  loop 
=  refrigerant  enthalpy  at  coil  exit  obtained  from  previous 
calculation  loop 
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ABSTRACT 


An  analysis  of  the  vapor  compression  cycle  and  the  main  components  of  an  air 
source  heat  pump  during  steady-state  operation  has  been  performed  with  emphasis 
on  fundamental  phenomena  taking  place  between  key  locations  in  the  refrigerant 
system.     The  basis  of  the  general  heat  pump  model  formulation  is  the  logic 
which  links  the  analytical  models  of  heat  pump  components  together  in  a  format 
reqiiiring  an  iterative  solution  of  refrigerant  pressure,  enthalpy  and  mass 
balances . 

The  modeling  effort  emphasis  was  on  the  local  thermodynamic  phenomena  which 
were  described  by  fundamental  heat  transfer  equations  and  equation  of  state 
relationships  among  material  properties.     In  the  compressor  model  several 
refrigerant  locations  were  identified  and  the  processes  taking  place  between 
these  locations  accounted  for  all  significant  heat  and  pressure  losses.  Eva- 
porator and  condenser  models  were  developed  on  a  tube-by-tube  basis  where  per- 
formance of  each  coil  tube  is  computed  separately  by  considering  the  cross-flow 
heat  transfer  with  the  external  air  stream  and  the  appropriate  heat  and  mass 
transfer  relationships.     A  capillary  tube  model  was  formulated  with  the  aid  of 
Fanno  flow  theory. 

The  developed  heat  pump  model  has  been  validated  by  checking  computer  results 
against  laboratory  test  data  for  full  and  part  load  operation  for  the  cooling/ 
dehumidif ying  mode  as  well  as  the  heating  mode  under  frosting .conditions . 

Key  Words:     air  conditioner;  capillary  tube;  coil;  compressor;  condenser; 

evaporator;  expansion  device;  heat  pump;  modeling;  vapor 
compression  cycle. 
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PREAMBLE 

This  report  is  essentially  the  dissertation  of  Piotr  Domanski  which  was 
submitted  to  the  Catholic  University  of  America  in  partial  fulfillment  of  his 
Ph.D.  requirements.    David  Didion  was  his  major  professor  and  dissertation 
director . 
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DISCLAIMER 


In  view  of  the  presently  accepted  practice  of  the  building  industry  in  the 
United  States  and  the  structure  of  the  computer  software  used  in  this  project, 
common  U.S.  units  of  measurement  have  been  used  in  this  report.     In  recognition 
of  the  United  States  as  a  signatory  to  the  General  Conference  of  Weights  and 
measures,  which  gave  official  status  to  the  SI  system  of  units  in  1960,  appro- 
priate conversion  factors  have  been  provided  in  the  table  below.     The  reader 
interested  in  making  further  use  of  the  coherent  system  of  SI  units  is  referred 
to:     NBS  SP330,  1972  Edition,   "The  International  System  of  Units,"  or  E380-72, 
ASTM  Metric  Practice  Guide  (American  National  Standard  2210.1). 


METRIC  CONVERSION  FACTORS 


Length 


1  inch  (in)  =  25.4  millimeters  (mm) 
1  foot  (ft)  =  0.3048  meter  (m) 


Area 


1  ft2  =  0.092903  m2 


Volume 


1  ft3  =  0.028317  m3 


Temperature 


F  =  9/5  C  +  32 


Temperature 
Interval 


1°F  =  5/9°C  or  K 


Mass 


1  pound  (lb)  =  0.453592  kilogram  (kg) 


Mass  Per  Unit 
Volume 


1  lb/ft3  =  16.0185  kg/m^ 


Energy 

Specific  Heat 
Gallon 


1  Btu  =  1.05506  kilojoules  (kJ) 

1  Btu/[(lb)(''F)]  =  4.1868  kJ/[(kg)(K)] 

1  gallon  =  0.0037854  m^ 
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1 .  INTRODUCTION 


Any  power  producing  cycle,  whose  processes  can  be  performed  in  the  reverse 
direction  and  in  the  reverse  order,  will  serve  as  a  refrigerating  cycle. 
The  Carnot  cycle  is  the  model  of  perfection  for  a  power  producing  cycle  as 
well  as  for  a  refrigeration  cycle.     However,  it  is  not  practical  to  realize 
this  perfect  thermodynamic  cycle  with  hardware  (see  figure  1).     First,  two 
compression  machines,  one  realizing  isentropic  compression,  the  other  realiz- 
ing isothermal  compression  would  be  required.     Second,  the  small  amount  of 
expansion  work  available  does  not  support  the  idea  of  incorporating  a  work 
producing  expander  into  the  system. 


CONDENSER 


i  W  =  WORK  TRANSFER 

jQg,  Q  =  HEAT  TRANSFER 


EMTROPY 


Figure  1 . 


Carnot  vapor  compression  cycle 


Thus  a  simple  isenthalic  throttling  expansion  device  and  only  one  type  of 
compressor  are  used  in  a  normal  refrigeration  system.     This  changes  the 
theoretical  cycle  as  shown  in  figure  2.     This  cycle  will  have  a  coefficient 
of  performance  smaller  than  that  for  the  Carnot  cycle. 

Figure  2  identifies  the  main  elements  of  a  system  working  on  the  vapor 
compression  cycle  principle.     They  are:     compressor,  condenser,  evaporator  and 
throttling  valve.     There  are  many  concepts  and  designs  of  heat  pump  components, 
however,  certain  types  are  predominant.     A  reciprocating  hermetic  compressor 
is  mostly  used  for  vapor  compression  in  small  systems.     Heat  exchangers  are 
usually  in  the  form  of  staggered  tubes  with  closely  packed  wavy  fins  (spine- 
fins  or  bristle-fins  are  sometimes  used).     There  are  basically  two  types  of 
expansion  devices  in  application  today.     Constant  flow  area  restrictors  are 
distinctly  prevailing  in  small  capacity  units  such  as  household  refrigerators, 
window-type  air  conditioners  and  central  residential  air  conditioners  and  heat 
pumps.     Larger,  more  expensive  residential  units  and  commercial  units  are  mostly 
equipped  with  a  variable  flow  area  device  called  a  thermostatic  expansion  valve 
(TXV) . 

The  focus  of  this  study  was  the  mathematical  modeling  and  prediction  of 
performance  of  the  vapor  compression  refrigeration  system  commercially  known 
as  a  residential  air  source  heat  pump.     The  type  of  heat  pump  considered  is 
the  one  most  commonly  available  with  a  reciprocating  hermetic  compressor,  fin 
tube  heat  exchangers,  and  a  constant  flow  area  expansion  device. 

The  present  report  discusses  and  analytically  describes  the  phenomena  taking 
place  in  the  heat  pump  system.     Based  on  this  analysis,  mathematical  models 
of  the  main  components  of  the  heat  pump  were  developed  and  linked  together  in 
a  general  heat  pump  simulation  model  HPSIM.     The  model  HPSIM  simulates  the 
performance  of  a  heat  pump  with  a  capillary  tube  at  given  outdoor  and  Indoor 
air  conditions  with  no  prefixed  assumption  on  refrigerant  state  at  any  point 
of  the  system. 
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Figure  2.     Theoretical  vapor  compression  cycle 
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2.     HEAT  PUMP  THERMODYNAMIC  CYCLE 


The  purpose  of  a  heat  pump  Is  to  pump  heat  from  a  low  temperature  reservoir  to 
a  high  temperature  reservoir.     From  the  thermodynamic  point  of  view,  the  basic 
heat  pump  components  needed  to  make  this  heat  transfer  feasible  are:     two  heat 
exchangers,  a  compressor,  and  an  expansion  device.     A  heat  pump  thermodynamic 
cycle  can  be  explained  considering  the  processes  that  the  refrigerant  undergoes 
in  these  four  components.     The  most  convenient  diagram  for  such  explanation  and 
performance  analysis  is  that  of  a  pressure  vs.  enthalpy  coordinate  system,  as 
shown  in  figure  3.     The  compressor  receives  low  pressure  and  temperature  refri- 
gerant at  state  1  and  compresses  it  to  a  high  pressure.     This  compression 
process  is  associated  with  an  increase  of  refrigerant  temperature.     At  state  2, 
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Figure  3o     Oversimplified  thermodynamic  cycle  of  a  heat  pump 
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the  high  pressure  and  high  temperature  vapor  enters  the  condenser.  The 
refrigerant  passing  through  the  condenser  rejects  heat  to  the  high  temperature 
reservoir  and  changes,  usually,  to  a  subcooled  liquid  at  state  3.     Then,  the 
refrigerant  flows  through  the  expansion  device  undergoing  a  drop  in  pressure 
and  temperature.     Finally,  the  low  pressure,  low  temperature,  and  low  quality 
refrigerant  at  state  4  enters  the  evaporator,  where  it  picks  up  heat  from  the 
low  temperature  reservoir,  reaching  a  superheated  (or  high  quality)  vapor 
state  1  at  the  evaporator  exit.     In  the  explanation  above,  the  low  and  high 
temperature  reservoirs  are  the  indoor  and  the  outdoor  environment,  when  the 
heat  pump  is  operating  in  the  cooling  mode. 

Besides  the  compressor,  two  heat  exchangers,  and  the  expansion  device,  there 
are,  for  practical  reasons,  many  other  components  in  an  actual  heat  pump  sys- 
tem.    Four  modeling  consideration,  the  most  important  are:     tubes  connecting 
basic  elements,  4-way  valve  enabling  refrigerant  flow  reversal  for  the  unit  to 
operate  in  the  heating  or  cooling  mode,  and  an  accumulator,  which  acts  as 
a  protective  device  for  the  compressor  by  storing  excess  refrigerant  during 
part  load  operation  and  preventing  liquid  refrigerant  from  entering  the  com- 
pressor.    All  the  above  mentioned  components  make  up  an  actual  vapor  compres- 
sion system.     Configuration  of  such  a  system  and  the  thermodynamic  cycle  are 
illustrated  in  figure  4  and  figure  5,  respectively.     Refrigerant  states  in 
figure  5  correspond  to  particular  locations  in  the  system  marked  in  figure  4, 
which  are: 

1  -  evaporator  exit,  suction  pipe  inlet 

2  -  suction  pipe  outlet,  low  pressure  4-way  valve  inlet 

3  -  low  pressure  4-way  valve  exit,  compressor  can  inlet 

4  -  inside  compressor  can 

5  -  compressor  cylinder  at  suction 

6  -  compressor  cylinder  at  discharge 

7  -  discharge  manifold 

8  -  compressor  can  exit,  high  pressure  4-way  valve  inlet 

9  -  high  pressure  4-way  valve  exit,  discharge  pipe  inlet 

10  -  discharge  pipe  outlet,  condenser  inlet 

11  -  condenser  outlet,  liquid  line  inlet 

12  -  liquid  line  outlet,  expansion  device  inlet 

13  -  evaporator  inlet 
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Note:     1.     Refrigerant  flow  direction  is  marked  as  for  cooling  operation 

2.     Numbers  5  and  6  situated  in  the  compressor  correspond  to 
refrigerant  state  before  and  after  compression  process 

Figure  4.     Schematic  of  a  heat  pump 
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Figure  5.     Thermodynamic  cycle  realized  by  a  heat  pump 
(For  number  location  refer  to  figure  4) 
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3.     OPERATION  OF  AN  AIR  SOURCE  HEAT  PUMP  WITH  A  CAPILLARY  TUBE 


Refrigerant  state  in  key  locations  of  a  heat  pump  and  heat  pump  performance 
are  functions  of  many  factors.     From  a  macroscopic  point  of  view  they  can  be 
listed  as  follows: 

-  individual  performance  of  heat  pump  components 

-  compatibility  of  components  in  the  system 

-  properties  of  refrigerant 

-  refrigerant  charge 

-  outdoor  and  indoor  air  conditions 
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Each  of  the  above  factors  have  an  independent  Impact  on  heat  pump  capacity  and 
efficiency. 

In  this  report,  a  properly  charged,  commercially  available  heat  pump  is 
discussed.     A  customer  receives  a  heat  pump  correctly  precharged  by  a  manu- 
facturer and  information  about  proper  charge  by  weight  of  refrigerant  is  sup- 
plied on  the  heat  pump  name  plates.     Simple  charts  are  also  available  to 
facilitate  checking  the  state  of  charge  of  unit  by  reading  compressor  suction 
and  head  pressures  at  given  operating  conditions. 

Manufacturers  establish  proper  charge  for  a  given  heat  pump  during  a  laboratory 
test.     It  is  usually  a  cooling  test  at  standard  operating  conditions,  which 
are : 

outdoor  air  dry  bulb  temperature  -  95°F 
indoor  air  dry  bulb  temperature  -  80°F 
indoor  air  wet  bulb  temperature    -  67 °F 

The  measure  of  the  refrigerant  charge  of  the  system  is  refrigerant  superheat 
at  the  compressor  can  inlet.     Superheat  of  about  15°F  at  this  location  is 
commonly  accepted  as  an  indication  of  a  proper  amount  of  refrigerant  in  the 
heat  pump. 

In  a  heat  pump  with  a  capillary  tube,  vapor  superheat  at  the  compressor  inlet 
will  vary  with  changing  operating  conditions.     As  a  matter  of  fact,  no  single 
refrigerant  parameter  stays  constant  when  indoor  or  outdoor  air  conditions  are 
altered  and  the  system  seeks  equilibrium  under  new  conditions.     To  understand 
the  phenomena  that  cause  refrigerant  state  changes,  consider  a  heat  pump  work- 
ing in  the  cooling  mode  at  constant  indoor  air  conditions  and  the  system 
response  to  a  step  increase  of  the  outdoor  temperature.     For  the  purpose  of  the 
analysis  it  is  assumed  that  the  overall  heat  transfer  coefficient  for  both  heat 
exchangers  does  not  vary.     It  is  also  assumed,  as  a  first  approximation,  that 
the  latent  heat  of  evaporation,  ifg,  is  constant  under  the  conditions  discussed. 

A  step  change  in  the  outdoor  air  temperature  first  affects  the  performance  of 
the  condenser.     Since  the  condenser  environment  temperature  increase,  the  mean 
temperature  difference  between  the  air  and  the  refrigerant  decreases.  Less 
heat  is  rejected  by  the  refrigerant  to  the  air  and  a  smaller  enthalpy  change 
is  realized  in  the  coil.     Consecutively,  higher  enthalpy,  less  subcooled  liq- 
uid refrigerant  enters  the  capillary  tube.     The  capillary  tube  is  sensitive  to 
the  amount  of  subcooling  and  with  less  subcooling  the  mass  flow  rate  through 
the  capillary  tube  decreases.     Since  the  compressor  capacity  stays  unchanged, 
pressure  builds  up  in  the  condenser.     Thus  point  11  on  the  p  -  h  diagram  (refer 
to  figure  5)  moves  in  the  direction  of  higher  enthalpy  and  higher  pressure. 

The  environment  of  the  evaporator  was  not  altered.     However,  refrigerant 
parameters  in  the  evaporator  change  in  response  to  the  change  of  refrigerant 
state  at  location  11  as  well  as  to  the  change  in  refrigerant  mass  flow  rate. 
For  pressure  drop  considerations,  as  a  first  approximation,  the  liquid  line, 
the  capillary  tube,  and  the  evaporator  are  viewed  as  one  tube,  in  which  a  given 
flow  experiences  a  certain  pressure  drop.     As  pressure  of  the  refrigerant  at 
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location  11  increases,  it  pulls  up  the  refrigerant  pressure  in  the  evaporator. 
Conflicting  phenomena  affect  the  change  of  refrigerant  enthalpy  at  point  1. 
Increased  enthalphy  at  point  11  and  reduced  refrigerant  mass  flow  rate  work 
against  the  impact  of  the  smaller  temperature  difference  between  the  air  and 
the  refrigerant  causing  refrigerant  enthalpy  to  slightly  decrease.     The  signi- 
ficance of  the  change  of  the  refrigerant  state  at  location  1  depends  on  the 
move  towards  smaller  vapor  specific  volume.     .Since  the  compressor,  as  a  first 
approximation,  is  a  constant  intake  volume  pump,  increase  in  gas  density 
results  in  a  higher  mass  flow  rate. 

The  indicated  phenomena  will  balance  themselves  further  until  the  refrigerant 
in  the  key  locations  of  the  system  acquires  thermodynamic  states  that  satisfy 
simultaneously  all  heat  pump  components.     Since  under  steady-state  operation, 
refrigerant  mass  flow  rate  through  the  compressor  and  the  capillary  tube  have 
to  match,  the  pressure  in  the  condenser  will  rise  to  an  appropriate  level. 
Again,  higher  pressure  in  the  condenser  implies  some  increase  of  pressure  (and 
saturation  temperature)  in  the  evaporator.     Smaller  temperature  difference 
between  the  indoor  air  and  the  refrigerant,  and  higher  refrigerant  mass  flow 
rate  result  in  a  smaller  refrigerant  enthalpy  increase  in  the  evaporator. 
Refrigerant  state  at  point  1  with  respect  to  point  11  is  determined  by  this 
enthalpy  change  and  appropriate  pressure  drop.     Finally,  when  steady-state 
conditions  are  obtained,  the  following  changes  in  refrigerant  thermodynamic 
states  can  be  noticed:     higher  refrigerant  pressure  in  the  condenser,  higher 
refrigerant  pressure  in  the  evaporator,  less  superheat  in  the  evaporator  exit, 
and  less  refrigerant  liquid  subcooling  at  the  condenser  exit.     A  change  of 
saturation  temperature  in  the  condenser  corresponds  approximately  to  a  change 
in  the  outdoor  air  temperature.     The  change  of  pressure  in  the  evaporator 
equals  to  a  fraction  (approximately  15  percent)  of  the  pressure  change  in  the 
condenser.     Though  the  refrigerant  is  circulated  in  the  system  at  a  higher 
mass  flow  rate,  capacity  of  the  heat  pump  is  decreased  due  to  smaller  refrig- 
erant enthalpy  change  in  the  evaporator.     Efficiency  also  drops  since,  in  addi- 
tion to  a  smaller  cooling  effect,  energy  input  for  the  compression  process 
increases  due  to  higher  refrigerant  mass  flow  rate  and  higher  compression  rate. 
A  decrease  of  the  outdoor  air  temperature  would  induce  the  opposite  trends. 

A  change  of  refrigerant  parameters  is  illustrated  in  figure  6,  where  results 
of  three  heat  pump  tests  are  plotted  on  the  pressure-enthalpy  diagram.  (Test 
results  are  those  of  system  2  discussed  later  in  section  6.1).     In  these  tests, 
the  indoor  air  conditions  were  held  constant  while  the  outdoor  air  dry  bulb 
temperature  was  changed  inducing  modification  of  the  realized  thermodynamic 
cycle. 

For  the  heating  mode  operation,  refrigerant  flow  is  reversed  by  reset  of  the 
4-way  valve.     Flow  direction  of  the  refrigerant  in  the  system  is  the  opposite 
to  that  in  the  cooling  mode,  with  the  exception  of  the  compressor.     The  indoor 
coil  becomes  a  condenser  while  the  outdoor  coil  serves  as  an  evaporator.  It 
is  worthwhile  to  note  the  effect  this  change  has  on  the  system.     Since  pres- 
sures in  heat  pump  heat  exchangers  are  functions  of  environment  temperatures, 
the  pressure  in  the  evaporator  is  now  much  lower  than  during  cooling  operations. 
Consequently,  the  density  of  the  refrigerant  at  the  evaporator  exit  is  smaller 
and  less  refrigerant  is  being  pumped  by  the  compressor.     In  order  to  maintain  a 
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Indoor  conditiono  (all  three  tests):     air  dry  bulb  temperature  -  80°F 

air  wet  bulb  temperature  -  67°F 

Outdoor  dry  bulb  temperature:    105 °F 

 95°F 

  32°F 

Note:     1.     Saturation  lines  have  been  generated  by  means  of  equations 
presented  in  appendix  A. 

2.     Refrigerant  states  in  a  heat  pump  were  obtained  from 

laboratory  tests  2.1,  2.2.,  and  2.3  discussed  in  section  6.1. 

Figure  6.     Simplified  thermodynamic  cycle  realized  by  a  heat  pump  in  the 
cooling  mode  at  different  outdoor  conditions 
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proper  pressure  in  the  condenser,  another,  properly  sized  capillary  tube  is  now 
used  in  the  system,  while  the  one  used  in  the  cooling  mode  is  bypassed. 

An  indoor  coil  is  usually  much  smaller  than  an  outdoor  coil.     Need  for  good 
dehuraidif ication  performance  is  the  primary  reason.     The  indoor  coil  working 
as  the  condenser  cannot  hold  as  much  refrigerant  as  the  outdoor  coil  during 
cooling  operation.     Since  the  condenser  holds  the  bulk  of  working  refrigerant 
in  the  system,  excess  refrigerant  has  to  be  stored  somewhere  during  heating 
operation.     For  storage  purposes  in  heat  pumps  equipped  with  a  capillary 
tube,  an  accumulator  is  provided. 

A  schematic  of  an  accumulator  is  given  in  figure  7.     The  accumulator  is  a  closed 
container  with  two  pipes.     The  longer,  bent  pipe  has  a  small  diameter  hole  on 
the  bottom  side,  and  is  connected  to  the  compressor.    The  other  pipe  connects 
the  tank  with  the  evaporator  exit. 

If  superheated  vapor  leaves  the  evaporator,  vapor  only  is  contained  in  the 
accumulator  and  no  special  function  is  fulfilled  by  the  accumulator  in  this 
case.     If  wet  vapor  enters  the  accumulator,  liquid  droplets  accumulate  on  the 
bottom    of  the  tank,  while  saturated  vapor  enters  the  pipe  leading  towards  the 
compressor.     Some  liquid  refrigerant  enters  the  vapor  stream  through  the  hole 
in  the  pipe,  driven  by  the  static  pressure  difference  in  the  liquid-vapor 
stream  interface.     It  should  be  noted  that  during  steady-state  operation  with 
liquid    in  the  accumulator,  qualities  of  the  vapor  entering  and  leaving  the 
accumulator  have  to  be  equal.     Since  the  small  diameter  hole  allows  for  a  small 
rate  of  liquid  flow  into  the  saturated  vapor  stream,  vapor  taken  in  by  the 
compressor  is  close  to  its  saturation  state.     Such  a  situation  at  the  compres- 
sor inlet  prevails  in  the  heating  mode  since,  with  a  possible  exception  to 
operation  at  high  outdoor  temperature,  there  is  always  liquid  refrigerant 
stored  in  the  accumulator. 

This  discussion  has  been  for  the  heat  pump  employing  a  constant  flow  area 
expansion  device.     Thermostatic  expansion  valves,  TXV's,  are  variable  opening 
orifice  devices.     The  amount  of  opening  is  varied  according  to  some  control 
function,  usually  constant  vapor  superheat  at  the  evaporator  exit.     In  heat 
pumps  employing  TXV's,  vapor  superheat  at  the  compressor  inlet  does  not  depend 
on  operating  conditions  but  rather  has  a  constant  preset  value. 
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Figure  7 .     Schematic  of  an  accumulator 
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4.     PREVIOUS  WORK  IN  STEADY-STATE  VAPOR  COMPRESSION  CYCLE  MODELING 


Modeling  of  the  vapor  compression  cycle  has  to  be  performed  by  modeling 
essential  processes  taking  place  in  the  system.     A  review  of  state-of-the-art 
steady-state  vapor  compression  cycle  modeling  will  be  conducted  at  two  levels. 
On  one  level,  presently  available  heat  pump  models  will  be  discussed.  On 
another  level,  the  discussion  will  deal  with  the  available  models  of  particular 
components  of  a  heat  pump.     The  latter  discussion  will  take  place  in  later 
sections . 

Although  some  equipment  manufacturers  have  highly  developed  proprietary  models, 
one  of  the  first  sophisticated  steady-state  models  of  the  vapor  compression 
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cycle  was  published  at  Massachusetts  Institute  of  Technology  [1].     The  MIT 
model  simulates  a  heat  pump  in  the  heating  mode  equipped  with  a  thermostatic 
expansion  valve,  TXV,  maintaining  a  constant  refrigerant  vapor  superheat  at  the 
evaporator  exit.     Input  data  required  to  the  program,  besides  heat  pump  data, 
are  indoor  air  conditions,  saturation  refrigerant  temperature  at  the  evaporator 
exit,  and  refrigerant  vapor  superheat  at  the  same  location.     Heat  pump  capacity, 
power  input,  and  outdoor  air  temperature  are  obtained  as  results  of  a  simula- 
tion run.     If  heat  pump  performance  data  at  given  operating  conditions  are 
required,  several  simulation  runs  have  to  be  performed  in  order  to  be  able  to 
interpolate  to  get  desired  results. 

The  Oak  Ridge  National  Laboratory  developed  a  heat  pump  simulation  model  [2] 
based  on  the  study  performed  at  MIT.     The  ORNL  program  includes  a  new  com- 
pressor model,  more  suitable  for  use  by  an  averaging  user.     In  addition  to 
a  thermostatic  expansion  device  available  in  the  MIT  version,  a  capillary  tube 
model  has  been  incorporated,  hence  a  heat  pump  equipped  with  either  type  of  an 
expansion  device  can  be  simulated  by  the  ORNL  program,     the  general  structure 
of  MIT  model  has  been  retained.     Input  data  to  the  program  include  refrigerant 
saturation  temperature  and  refrigerant  vapor  superheat  at  the  evaporator  outlet 
while  as  in  the  MIT  model,  outdoor  air  temperature  is  one  of  the  outputs  of  the 
programl./ .     A  limitation  appears  in  the  simulation  of  a  system  employing  a  con- 
stant flow  area  expansion  device.     For  a  heat  pump  with  a  capillary  tube,  super- 
heat is  usually  known  only  at  one  set  of  standard  operating  conditions  (function 
of  refrigerant  charge  only).     Then,  when  a  change  in  operating  conditions  occurs, 
the  superheat  will  also  change.     The  ORNL  model  does  not  provide  a  means  for 
calculating  the  change  in  superheat;  it  is  required  as  input  data,  which  can 
be  approximated  at  best.     An  additional  limitation  also  results  from  the 
capillary  tube  model.     The  model  depends  on  curve  fitting  of  capillary  tube 
capacity  charts  presented  in  ASHRAE  Guide  and  Data  Book  [4].     These  charts 
cover  performance  of  a  long  capillary  tube  with  subcooled  liquid  entering. 
Thus ,  the  model  cannot  handle  either  a  capillary  tube  when  there  is  incomplete 
condensation  in  the  condenser  (that  does  happen,  particularly  during  heating 
operation),  or  very  short  constant  bore  flow  restrictors.     A  capillary  tube 
model  that  could  handle  all  cases  described  above  does  not  exist  in  literature 
to  the  author's  knowledge.     The  ORNL  model  was  validated  in  the  heating  mode 
only. 

A  mathematical  model  of  a  heat  pump  was  also  developed  at  the  National  Bureau 
of  Standards  [3].     The  model  simulates  the  performance  of  a  heat  pump  employing 
a  TXV  as  an  expansion  device.     The  general  logic  of  the  model  is  similar  to  those 
of  the  MIT  and  ORNL  models.     Heat  pump  data,  indoor  air  conditions,  refrigerant 
saturation  temperature,  and  refrigerant  vapor  superheat  at  the  evaporator  exit 
are  required  as  input.     Output  consists  of  performance  information  (capacity. 


±!  At  the  stage  of  completion  at  this  report  documentation  on  a  new  heat  pump 
model  was  published  by  Oak  Ridge  National  Laboratory  (December  1981)  [55]. 
This  latest  model,  representing  an  upgraded  version  of  the  model  described 
in  [2],  allows  the  user  to  specify  heat  pump  operating  conditions. 
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power  input,  efficiency)  and  the  outdoor  temperature.     The  NBS  model  was 
validated  in  the  heating  mode  only. 

The  above  review  outlines  the  available,  general  models  of  a  heat  pump,  their 
applicability  and  limitations.     Considering  the  formulation  of  a  model  of  a 
heat  pump  with  a  capillary  tube,  the  following  conclusions  can  be  made: 

1.  The  only  program  available  in  the  literature  for  simulating  a  heat  pump 
with  a  capillary  tube  is  provided  by  ORNL.     The  program  has  limitations 
which  result  from  the  logic  of  the  program  and  from  the  capillary  tube 
model  itself.     As  input  data,  refrigerant  vapor  superheat  at  the  evaporator 
exit  is  required,  and  the  program  is  unable  to  simulate  the  heat  pump  when 
a  two-phase  flow  condition  exists  at  the  capillary  tube  inlet.     Also,  the 
model  cannot  handle  a  short  capillary  tube. 

2.  The  logic  in  the  existing  models  is  adequate  for  simulation  of  a  heat  pump 
with  TXV,  however,  it  is  insufficient  for  modeling  a  system  employing  a 
capillary  tube.     A  new  logic  allowing  for  calculation  of  the  vapor 
superheat  should  be  formulated. 

3.  There  is  no  capillary  tube  model  available  in  the  literature  which  can 
simulate  its  performance  regardless  of  refrigerant  state  (two-phase  or 
subcooled  liquid)  at  the  tube  inlet,  and  which  can  simulate  tubes  of  a 
short  length.     Such  a  model  should  be  formulated. 

4.  Models  of  a  compressor,  an  evaporator,  and  a  condenser  are  available  in 
the  literature.     They  may  be  used  in  original  or  modified  form  in  any  new 
simulation  heat  pump  model. 

5.  None  of  the  available  models  accepts  operating  conditions  as  input  data. 

A  model  providing  heat  pump  performance  data  at  given  operating  conditions 
is  most  desirable  and  appropriate  logic  should  be  developed. 

6.  All  discussed  models  have  been  verified  in  the  heating  mode  only. 
Verification  of  a  model  in  heating  and  cooling  operation  would  be  desirable. 

These  conclusions  led  to  the  development  of  the  new  model  described  in  this 
report.     It  was  felt  that  the  model  should  simulate  an  air-to-air  heat  pump 
equipped  with  a  constant  flow  area  expansion  device  with  no  prefixed  refrigerant 
state  in  any  system  location  required  in  a  set  of  input  data. 
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5.     MODELING  A  HEAT  PUMP  WITH  A  CONSTANT  FLOW  AREA  EXPAI^SION  DEVICE 


5.1    MODELING  OF  THE  SYSTEM 

5.1.1    Structure  of  Available  Heat  Pump  Models 

During  heat  pump  operation,  refrigerant  parameters  at  any  location  are 
established  at  a  level  which  satisfies  the  respective  components  of  the  system 
There  is  a  one-to-one  relationship  between  working  medium  parameters  and  opera 
ing  conditions,  i.e.,  for  given  outdoor  and  indoor  air  conditions,  there  is  no 
alternative  but  just  one  set  of  refrigerant  parameters  at  any  location  within 


21 


the  system  which  constitute  steady  state  operation.     This  unique  set  of 
refrigerant  parameters  has  to  be  sought  by  a  heat  pump  simulation  program  in 
an  iterative  process. 

In  order  to  set  up  a  heat  pump  system  logic,  balances  taking  place  during 
steady  state  operation  have  to  be  recognized.     From  the  fact  that  a  thermody- 
namic cycle  is  a  closed  loop,  an  analysis  of  a  heat  pump  cycle  on  a  pressure- 
enthalpy  diagram  allows  statement  of  two  balances: 

-  Enthalpy  Balance 

This  balance  implies  that  energy  input  to  the  system,  i.e.,  to  the 
evaporator,  the  suction  line  and  the  compressor,  has  to  be  equal  to 
energy  rejected  by  the  system,  i.e.,  by  the  high  pressure  line  and 
the  condenser. 

-  Pressure  Balance 

This  balance  implies  that  increase  of  refrigerant  pressure  during  the 
compression  process  has  to  be  equal  to  the  total  pressure  drop  during 
other  processes  forming  the  cycle.     As  pressure  drop  and  mass  flow  rate 
are  strictly  related,  this  balance  may  be  restated  that  pressure  drop 
through  each  component  has  to  be  adequate  so  mass  flow  rate  in  each 
component  is  the  same. 

These  two  balances  are  sufficient  to  set  up  logic  for  a  program  simulating  the 
performance  of  a  heat  pump  equipped  with  a  thermostatic  expansion  valve.  This 
was  done  for  all  three  models  reviewed  [1],  [2],   [3].     The  logic  of  these  models 
is  very  similar,  so  only  one,  the  N.B.S.  model  logic  will  be  described. 

For  the  heating  mode  simulation,  besides  heat  pump  data,  indoor  air  conditions, 
refrigerant  saturation  temperature,  and  refrigerant  vapor  superheat  at  the 
evaporator  exit  are  required.     These  data  (including  refrigerant  state  at  the 
evaporator  exit)  are  fixed  for  the  run.     A  set  of  input  data  also  includes  com- 
pressor head  pressure  and  outdoor  air  temperature,  which  are  input  as  initial 
estimates. 

The  compressor  is  simulated  first  to  estimate  refrigerant  mass  flow  rate.  Then 
an  enthalpy  balance  is  sought  simulating  the  condenser  and  the  evaporator.  If 
refrigerant  enthalpy  at  the  condenser  outlet  and  the  evaporator  inlet  do  not 
match,  outdoor  air  temperature  is  altered  and  the  evaporator  is  simulated  until 
the  balance  is  reached.     It  is  assumed  reasonably  well,  that  a  change  in  outdoor 
temperature  has  negligible  effect  on  compressor  performance.     Next,  the  thermo- 
static expansion  valve  is  simulated.     If  refrigerant  mass  flow  rate  through  the 
TXV  is  not  equal  to  that  through  the  compressor,  the  compressor  discharge 
pressure  is  changed  and  the  whole  procedure  is  repeated  again. 

The  above  logic  is  relatively  simple,  since  all  iterative  processes  take  place 
with  a  fixed  refrigerant  state  at  the  evaporator  exit.     The  evaporator  environ- 
ment temperature  is  calculated  and  is  part  of  the  simulation  run  output.  If 
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performance  of  the  unit  at  particular  outdoor  conditions  is  needed,  a  few 
simulation  runs  have  to  be  performed  in  order  to  allow  for  interpolation.  Of 
course,  this  results  in  increased  computing  time  and  imposes  a  need  to  reduce 
output  data.     A  similar  logic  can  be  used  to  simulate  the  cooling  mode,  however, 
with  similar  inefficiencies. 

5.1.2    New  System  Logic 

Two  balances,  an  enthalpy  balance  and  a  pressure  balance,  are  sufficient  to  set 
up  the  logic  for  a  program  simulating  steady  state  performance  of  a  heat  pump 
with  a  thermostatic  expansion  valve  because  this  expansion  device  maintains  a 
constant  refrigerant  vapor  superheat  at  evaporator  outlet.     A  constant  flow 
area  expansion  device  does  not  provide  any  fixed  control  of  refrigerant  super- 
heat or  any  other  flow  parameter.     It  can  be  generally  stated,  that  for  heat 
pumps  with  capillary  tubes  during  cooling  mode  operation,  the  vapor  superheat 
at  the  compressor  can  inlet  increases  with  rising  indoor  temperature  and  humi- 
dity, and  with  decrease  of  outdoor  temperature.     Logically  then,  the  superheat 
decreases  with  respective  opposite  trends.     During  cooling  operation,  superheat 
may  vary  from  0  to  40°F.     During  heating  mode  operation,  very  low  or  no 
superheat  is  observed  at  the  compressor  inlet. 

In  as  much  as  refrigerant  superheat  is  a  new  unknown,  an  additional  balance  has 
to  be  used  to  establish  superheat  values  for  each  operating  condition.     It  is 
proposed,  that  superheat  be  evaluated  by  means  of  conservation  of  mass,  as 
expressed  by  the  equation: 

2-  /  pdV  =  0  (1) 
Dt  V 

where         p  =  density 

V  =  system  volume 

Equation  (1)  implies,  that  the  amount  of  working  medium  in  the  system  is  the 
same  at  all  times,  for  any  operating  conditions. 

The  logic  is  based  on  three  balances,  enthalpy,  pressure,  and  refrigerant  mass 
as  explained  below.     It  provides  means  for  calculation  of  vapor  superheat  and 
allows  for  heat  pump  simulation  at  implied  operating  conditions.  Explanation 
of  the  logic,  for  sake  of  clarity,  is  done  considering  only  four  main  components 
of  a  heat  pump,  i.e.,  a  compressor,  a  condenser,  a  capillary  tube,  and  an 
evaporator. 

Required  input  consists  of  outdoor  and  indoor  conditions,  and  heat  pump  data. 
The  simulation  process  starts  with  guessed  refrigerant  pressure  and  vapor 
superheat  at  the  compressor  inlet,  and  the  compressor  discharge  pressure. 
Using  these  data,  compressor  performance  is  simulated  yielding  refrigerant  mass 
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flow  rate.     Next,  the  condenser  and  the  capillary  tube  are  simulat  edl/  and  a 
mass  flow  balance  is  sought  by  comparing  refrigerant  mass  flow  rate  through  the 
compressor  and  capillary  tube.     If  compared  mass  flows  are  not  equal,  simula- 
tion of  the  compressor,  the  condenser,  and  the  capillary  tube  is  redone  with  an 
unchanged  refrigerant  state  at  the  compressor  can  inlet  and  a  modified  guess  of 
compressor  discharge  pressure.     Increasing  this  pressure  reduces  refrigerant 
mass  flow  rate  through  the  compressor.     Then  the  condenser  works  with  this 
smaller  mass  flow  rate  and  at  a  higher  saturation  temperature.  Consequently, 
refrigerant  reaching  the  expansion  device  inlet  is  at  higher  pressure  and  has 
more  subcooling.     Both  factors  promote  greater  mass  flow  rate  through  the  capil- 
lary tube.     Thus,  increase  in  discharge  pressure  has  a  clear  and  opposite  impact 
on  mass  flow  rates  through  the  compressor  and  through  the  capillary  tube  and 
the  appropriate  discharge  pressure  can  be  found     >r  which  mass  flow  balance 
exists . 

Once  mass  flow  balance  is  reached,  simulation  of  the  evaporator  is  performed 
with  the  known  refrigerant  state  at  the  evaporator  exit  and  the  refrigerant 
mass  flow  rate.     Since  the  thermodynamic  process  in  a  capillary  tube  may  be 
viewed  as  adiabatic,  refrigerant  enthalpy  at  the  evaporator  inlet  should  be 
equal  to  the  enthalpy  at  the  condenser  outlet.     If  these  enthalpies  are  not 
equal  (enthalpy  balance  is  not  reached),  a  new  calculation  starts  from  the 
beginning  with  a  modified  guess  on  the  refrigerant  pressure  at  the  compressor 
can  inlet.     From  the  condenser  operation  point  of  view,  a  change  in  compressor 
suction  presure  induces  a  change  in  refrigerant  mass  flow  rate  and  modification 
of  the  condenser  saturation  temperature  caused  by  the  system  flow  balance 
search.     These  two  changes  have  opposite  effects  on  refrigerant  enthalpy  at  the 
condenser  exit,  leaving  it  only  slightly  altered.     On  the  other  hand,  the  same 
change  of  compressor  suction  pressure  has  a  strong  effect  on  refrigerant 
enthalpy  at  the  evaporator  inlet  by  a  change  in  the  evaporator  saturation  tem- 
perature and  refrigerant  mass  flow  rate,  both  working  to  change  the  enthalpy  in 
the  same  direction.     Thus  an  appropriate  suction  pressure  at  the  compressor 
inlet  can  always  be  found  for  which  an  enthalpy  balance  exists. 

Once  enthalpy  and  pressure  balances  are  established,  two  out  of  three 
refrigerant  parameters  guessed  at  the  beginning  are  determined.     However,  these 
two  parameters,  compressor  suction  and  discharge  pressure,  are  only  as  good  as 
the  third  guessed  parameter,  the  refrigerant  superheat  at  the  compressor  inlet, 
which  has  not  yet  been  verified.     For  this  parameter  verification,  the  refrige- 
rant mass  balance  in  the  system  is  used.     Specifically,  a  refrigerant  mass 
inventory  is  made.     It  is  based  on  refrigerant  states  in  the  system  which  are 
found  to  satisfy  enthalpy  and  pressure  balances  with  asumed  vapor  superheat  at 
the  compressor  can  inlet.     Then  the  amount  of  refrigerant  obtained  from  mass 
inventory  calculations  is  compared  to  the  actual  refrigerant  mass  input  (known 


tl   Simulation  of  a  capillary  tube  can  be  performed  at  this  time  even  without 
precise  knowledge  of  the  evaporator  inlet  pressure.     For  a  heat  pump  pres- 
sure operation  range,  capillary  tube  capacity  is  a  weak  function,  if  not 
independent  of  evaporator  pressure.     For  the  first  calculation  pressure  at 
the  evaporator  inlet  may  be  assumed  to  be  the  same  as  at  the  evaporator 
outlet.     Refer  to  section  5.3. 
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as  design  parameter) .     If  the  amount  of  refrigerant  calculated  is  smaller  than 
refrigerant  input  into  the  heat  pump,  the  superheat  guess  has  to  be  decreased 
and  all  calculations  have  to  be  repeated  from  the  beginning. 

The  logic  described  above  is  presented  graphically  in  figure  8.     Details  about 
the  above  program  flow  chart  can  be  found  in  appendix  H.     Refrigerant  mass 
inventory  is  discussed  in  section  5.6. 

5.2    MODELING  OF  A  RECIPROCATING,  HERMETIC  COMPRESSOR 

5.2.1     Compressor  Operation 

The  compressor  is  mechanically  the  most  complex  component  of  a  heat  pump.  A 
reciprocating,  hermetic  compressor  is  most  commonly  used  in  heat  pump  systems. 
A  schematic  of  this  type  of  compressor  is  shown  in  figure  9.  The  compressor 
consists  of  a  shell  containing  an  electric  motor,  a  cylinder/piston  assembly 
with  valves  and  manifolds  and  pipes.  The  electric  motor  is  coupled  to  the 
compressor  eccentric.  Lubricating  oil  is  collected  on  the  bottom  of  the  can 
and  is  in  free  contact  with  the  refrigerant. 

Flow  of  the  refrigerant  is  from  location  3  to  location  8,  as  marked  in  figure  9. 
Low  pressure  refrigerant  enters  the  compressor  can  and  is  directed  towards  the 
electric  motor  to  cool  motor  windings.     Enthalpy  of  the  refrigerant  changes  due 
to  this  and  other  heat  transfer  with  other  surrounding  surfaces;  namely,  the 
compressor  can,  manifolds,  the  discharge  pipe  and  the  cylinder  body.     Then  the 
refrigerant  at  state  4  undergoes  the  cylinder  process,  which  in  case  of  posi- 
tive displacement  compressors  consists  of:     expansion  through  the  suction  valve 
and  mixing  the  residual  gas,  compression,  expansion  through  the  discharge  valve, 
and  re-expansion  of  residual  gas.     Compressed,  high  temperature  and  high  pres- 
sure refrigerant  vapor  leaves  the  compressor  through  the  discharge  manifold  and 
discharge  line  after  giving  up  some  heat  to  the  low  pressure  refrigerant  at 
state  4. 

Numbers  marked  in  figure  9  identify  key  compressor  locations  and  corresponding 
refrigerant  states.     Even  during  steady  state  operation  it  is  difficult  to  talk 
about  steady  parameters  in  the  compressor  can,  where,  pressure  changes  in  a 
manner  of  pulses.     This  pulsation  is  initiated  by  the  periodic  compression  and 
suction  processes  and  valve  behavior.     Because  a  certain  pressure  difference  is 
needed  for  valve  opening,  pressure  in  the  cylinder  at  the  end  of  compression  is 
considerably  above  discharge  manifold  pressure.     Similarly,  pressure  during 
intake  falls  below  suction  manifold  pressure.     Particularly  strong  pressure 
peaks  are  observed  in  case  of  spring  equipped  valves.     These  pressure  differ- 
ences vary  from  compressor  to  compressor  because  of  valve  and  manifold  design, 
but  also  vary  for  the  same  compressor  with  compressor  speed,  compression 
pressure  ratio,  refrigerant,  and  mass  flow  rate. 

Since  processes  taking  place  in  a  compressor  are  dynamic  in  nature,  the  most 
correct  way  of  modeling  a  compressor  would  require  dynamic  simulation  of  valve 
motion,  cylinder-manifold  pressure  interaction,  and  cylinder  heat  transfer. 
Such  models  exist  and  are  described  in  [5],   [6],  and  [7].     However,  these  models 
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Overall  logic  of  a  program  simulating  a  heat  pump  with 
constant  flow  area  expansion  device 
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usually  require  some  very  detailed  experimental  data  or  design  information,  not 
necessarily  readily  available  for  the  prospective  user  of  a  heat  pump  simulation 
program.     That  is  why  a  relatively  simpler,  but  certainly  still  adequate 
compressor  model  from  [3]  was  adopted  for  this  study. 

5.2.2    Theory  and  Governing  Relations 

Several  assumptions  are  made  for  model  formulation.     The  basic  assumption  is 
that  the  highly  dynamic  process  in  the  compressor  results  in  steady  gas  param- 
eters through  the  can.     Pressure  and  temperature  in  the  cylinder  and  in  suction 
and  discharge  manifolds  are  assumed  to  have  constant  values.     Refrigerant  is 
considered  to  have  uniform  thermodynamic  properties  throughout  the  space 
assigned  to  the  particular  location.     Refrigerant  flow  through  the  compressor 
is  assumed  to  be  one-dimensional,  so  one-dimensional  steady  flow  equations  can 
be  used.     Though  the  model  does  not  simulate  valve  dynamic  behavior,  it  allows 
for  valve  representation  by  a  steady  difference  between  manifold  and  cylinder 
pressures. 

In  order  to  derive  equations  governing  compressor  balances,  the  thermodynamic 
irreversibilities  taking  place  in  the  hermetically  sealed  compressor  have  to  be 
identified.     These  losses  can  be  put  into  four  categories: 

1.  Those  associated  with  incomplete  conversion  of  electric  energy  into  a 
mechanical  energy  available  tor  vapor  compression. 

2.  Those  associated  with  non-isentropic  conversion  process  inside  the 
cylinder. 

3.  Those  due  to  heat  transfer  at  different  locations. 

4.  Those  due  to  pressure  drop  at  difference  locations. 

All  of  these  losses  contribute  to  the  overall  compressor  efficiency  on  the 
order  of  60  percent  (first  law  basis)  and  have  to  be  considered  in  the  modeling 
effort . 

Conversion  of  Electrical  to  Mechanical  Energy 

Electric  energy  is  supplied  to  an  electric  motor  to  be  converted  into  mechanical 
work.     This  conversion  has  losses  due  to  windage,  friction,  winding  resistance, 
and  hysteresis,  which  are  accounted  for  by  a  motor  efficiency  rie*     By  definition 


(2) 


where 


E  =  electric  power  input 
Wg  =  mechanical  power  output 
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Electric  motor  efficiency,  ne>  depends  on  load  and  is  customarily  given  as  a 
function  of  the  fraction  of  actual  mechanical  power  output  to  the  maximum  power 
output.     In  heat  pumps  with  usual  power  requirements  below  5  hp,  single  phase 
electric  motors  are  mostly  used.     They  are  permanent  split-capacitor  or 
capacitor-start  capacitor-run  types.     A  typical  efficiency  versus  load  curves 
for  a  permanent  split-capacitor  2  pole  electric  motor  is  shown  in  figure  10 
[14]. 

The  electric  motor  is  coupled  to  the  compressor  eccentric.     Moving  compressor 
parts  experience  friction  and  a  power  loss  accounted  for  by  the  compressor 
mechanical  efficiency  nm*     By  definition 

e 

where  =  mechanical  power  available  for  compression 

A  common  value  of  compressor  mechanical  efficiency  falls  in  range  of  0.95  to 
0.98. 

Compression  Process 

A  simplified  Indicator  diagram  Is  shown  in  figure  11.     It  shows  assumed  constant 
suction  and  discharge  pressures  in  a  cylinder.     In  modeling  the  cylinder  process, 
both  compression  and  re-expansion  processes  are  assumed  to  be  polytropic  with 
the  same  polytropic  Index,  n,  following  the  equation: 

P  •  v'^  =  const  (4) 

where         n  =  polytropic  index 
P  =  pressure 

V  =  refrigerant  vapor  specific  volume 

Since  constant  pressure  and  temperature  are  assumed  during  the  discharge  process, 
this  Implies  no  change  in  specific  volume  between  points  "c"  and  "d"  (figure 
11);  consequently,  the  compression  and  re-expansion  polytropic  curves  will  coin- 
cide.    This  further  means,  that  the  net  work  required  for  compression  of  the 
residual  gas  is  zero. 

The  refrigerant  enthalpy  Increase  during  polytropic  compression,  I5  -  15  (refer 
to  figure  5),  can  be  evaluated  by  the  equation  derived  from  the  expressions  for 
Isentroplc  and  polytropic  work  of  compression  at  the  same  compression  ratio. 
Equating  these  expressions  results  in  the  equation: 


^    fe)^'^'  -1 


(5) 


Y  -  1  .P 


A]  -  1 

P5 
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Figure  10.     Typical  efficiency  versus  load  curve  for  a  permanent 
split-capacitor  2  pole  electric  motor  [4] 
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Figure  11.     Simplified  indicator  diagram  for  a  reciprocating  compressor 
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where        I5  -  refrigerant  enthalpy  before  compression 
I5  =  refrigerant  enthalpy  after  compression 
Igg  =  refrigerant  enthalpy  after  Isentroplc  compression,  defined  by 
pressure  P5  and  entropy  S5 
n  =  polytroplc  Index 
P5  =  suction  pressure 
P5  =  discharge  pressure 
Y  =  Isentroplc  Index 


The  Isentroplc  Index,  Y>  and  the  polytroplc  Index  n,  are  related  by  the 
polytroplc  efficiency  of  the  compressor: 


Y  -1 

P      n  -  1 
n 


(6) 


The  Isentroplc  Index,  Yj  equals  the  ratio  of  specific  heats,  as  In  the 
following  formula: 

Y  -  ^  (7) 
Cv 

where        Cp  -  specific  heat  at  constant  pressure 
Cv  =  specific  heat  at  constant  volume 

The  specific  heats  of  ideal  gases  are  constant,  thus  the  Isentroplc  index  of  an 
ideal  gas  is  also  constant.     For  real  gases  as  refrigerant  vapors,  specific 
heats  vary  along  the  compression  path.     To  accommodate  this  fact,  the  Isentroplc 
index  can  be  evaluated  by  taking  the  average  of  the  respective  specific  heat 
ratios  at  points  5  and  6s. 

The  refrigerant  enthalpy  Increase  during  polytroplc  compression  could  also  be 
calculated  using  isentroplc  efficiency  and  enthalpy  Increase  during  isentroplc 
compression  at  the  same  compression  ratio.     However,  isentroplc  efficiency  is 
sensitive  to  the  compression  pressure  ratio  while  polytroplc  efficiency  is  more 
consistent  from  one  application  to  another  and  provides  more  consistent  repre- 
sentation of  average  compressor  performance  [8].     The  imperfection  of  using 
isentroplc  efficiency  may  be  traced  to  the  general  thermodynamic  relation: 


=  T  (8) 

\  3s  /  p 

where         1  =  enthalpy 
P  =  pressure 
s  =  entropy 
T  =  temperature 


requiring  that  pressure  lines  diverge  on  a  Molller  chart. 
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The  refrigerant  mass  flow  rate  pumped  by  a  compressor  can  be  calculated  by  the 
following  formula: 


m^  =  60  »  RPM  '  Vs 

where        rnj-  =  refrigerant  mass  flow  rate 

RPM  =  compressor  speed  (number  of  revolutions  per  minute) 
Vs  =  compressor  swept  volume  per  revolution 

V5  =  refrigerant  specific  volume  at  piston  before  compression 
Tiy  =  volumetric  efficiency 


The  RPM  of  a  compressor  is  equal  to  that  of  an  electric  motor  and  is  a  function 
of  load  on  the  motor.     A  typical  speed  (RPM)  versus  load  curve  for  a  permanent 
split-capacitor  2-pole  elctric  motor  is  shown  in  figure  12. 

With  assumptions  stated  so  far  for  the  compression  and  re-expansion  processes, 
the  following  formula  for  volumetric  efficiency  of  a  compressor  can  be  derived: 


1 


(10) 


where        C^.  =  correction  factor  for  leakage  from  piston,  valves  and  for 
throttling  effect,  assumed  to  be  0.96,  [9] 
Cg  =  clearance  volume,  in  fraction  of  stroke  volume 


Heat  Transfer  Relations 


Electrical  energy  supplied  to  an  electric  motor  is  in  part  transferred  to  the 
refrigerant  and  in  part  is  dissipated  to  the  compressor  environment.  This 
compressor  heat  balance  can  be  expressed  by  the  following  equation: 


E  +  mr  (i3  -  is)  "  Qc,A  =  0  (H) 

where  E  =  electric  power  input 

13,13  =  refrigerant  enthalpy  at  respective  locations  (refer  to  figure  9) 
Qq^P^  =  rate  of  heat  rejected  to  ambient  air 

Reference  is  made  to  figure  9,  where  key  locations  of  refrigerant  in  the  hermetic 
compressor  are  marked.     In  order  to  solve  the  compressor  heat  balance  represented 
by  equation  (11),  the  following  heat  transfer  losses  are  considered  in  the 
compressor  model: 

1.  Heat  transfer  between  the  compressor  can  and  ambient  air,  Qc,a 

2.  Heat  transfer  between  the  refrigerant  in  the  compressor  can  and  the 
compressor  can,  Q4^c 
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Figure  12.    Typical  speed  (RPM)  versus  load  curve  for  a  permanent 
split-capacitor  2  pole  electrical  motor  [4]. 
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3.  Heat  transfer  between  the  inlet  refrigerant  and  the  suction  toanifold 
and  valve,  Q4^5 

4.  Heat  transfer  between  the  discharge  refrigerant  and  the  discharge 
valve  and  manifold,  Q6^7 

5.  Heat  transfer  between  the  refrigerant  in  the  compressor  can  and  the 
refrigerant  in  the  discharge  line,  Qy^s 

According  to  general  heat  transfer  equation  of  P^clet: 

Q  =  U  •  Ah  •  AT  (12) 

The  heat  transfer  rate,  Q,  for  a  given  temperature  difference,  AT,  can  be 
calculated  if  the  overall  coefficient  of  heat  transfer,  U,  and  heat  transfer 
area,  A|^,  are  known.     The  heat  transfer  between  the  compressor  can  and  ambient 
air  (item  1)  is  governed  by  free  convection  (radiation  is  negligible).     A  simple 
correlation  is  available  to  estimate  the  free  convection  heat  transfer 
coefficient  in  air  [10],  i.e., 

h  a  AtO-333  (13) 

where         h  =  convection  heat  transfer  coefficient 

AT  =  temperature  difference  between  air  and  heat  transfer 
surface 

Hence,  the  heat  transfer  rate  between  ambient  air  and  the  compressor  shell, 
Qc,A>  ^®  calculated  by  the  equation: 

QC,A  =  CQc,A(Te  -  T^)^-^^^  d^) 

where  CQc  a  ~  heat  transfer  parameter 
T^  =  ambient  air  temperature 
Tj,  =  compressor  shell  temperature 

Heat  transfer  inside  the  compressor  shell  between  the  shell  and  refrigerant 
vapor  (item  2)  is  governed  by  forced  convection,  for  which  the  non-dimensional 
heat  transfer  parametric  expression  in  terms  of  Nusselt,  Reynolds  and  Prandtl 
numbers  is  in  the  form: 

Nu  «  ReO*8  •  PrO*333  (15) 

where        jju  =  ^  *  ^         Nusselt  number 
k 

Pr  =  y  *  Prandlt  number 

k 

Re  =  ^  '  Reynolds  number 
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G  =  refrigerant  mass  flux 
A 

L  =  characteristic  dimension 

k  =  refrigerant  thermal  conductivity 

y  =  refrigerant  dynamic  viscosity 

Using  equation  (15)  the  forced  convection  heat  transfer  coefficient,  h,  can  be 
expressed  as: 

h  =  C  .  m^O'8  ,  ^0.667  .  .  ^  -  0.467 

where  C  =  constant  of  proportionality,  a  fuction  of  wetted  surface 

geometry 

Combining  equations  (12)  and  (16)  yields  the  following  expression  for  heat 
transfer  rate  between  the  compressor  can  and  refrigerant: 

n  -PA  ™  0»8  1,  0.667  ^„  0.333  ,  -  0.467  _  t  ^ 
%,C  ~  ^  '  \  '  •  "^4  •  ^P4  '  H  ^^4  T^) 

or 

A        —  rr\  0.8      1    0.667  0.333  —  0.467  /rp    _  t  \  /'^7^ 

^4,0  ~  '-'^4,0  *  ™r        *  H  •  ^P4  *  1^4  ^c* 

where  CQ4^c  ~  heat  transfer  parameter 

Obviously,  heat  transfer  rates        ^  and  Q4  q  are  equal.     Derivations  of  equations 
for  the  heat  transfer  rates  between  the  iniet  refrigerant  and  the  suction  mani- 
fold and  valve,  Q4  5,  and  between  the  discharge  refrigerant  and  the  discharge 
manifold  valve,  Q5  7,  is  similar  to  the  derivation  just  performed  for  Q4  q,  since 
in  these  two  cass  heat  transfer  is  also  by  forced  convection.  Resulting 
expressions  for  the  heat  transfer  rates  are: 

r.        -  rn  n,  0«8      .       0.666  0.333  -  0,467  _  t  ^ 

(18) 

n        -  rn  ™  0'8      1.      0.667      p„      0.333  -  0.467  /rp    _  t  ^ 

^6,7  -  CQ6,7  •  °^r        •  ^,7  *  ^P6,7  '  %  ,7  ^^7 

Heat  transfer  between  refrigerant  in  the  compressor  can  and  refrigerant  in  the 
discharge  line  (item  5)  is  modeled  as  forced  convection  heat  transfer  between 
the  fluids  separated  by  a  barrier  non-resistant  to  heat  flow.     Assuming  that 
the  temperature  of  refrigerant  in  the  shell  does  not  change  (as  a  result  of 
other  heat  transfers  in  the  can  and  mixing)  and  applying  the  logarithmic  mean 
temperature  difference,  the  following  expression  for  heat  transfer  rate  Q7  3 
can  be  derived: 

Q7,8  -  Cq7.8  •  (T7  -  T8)/l„Il^  (20) 
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where        al  =  (Cp^  •  Cp^  ^g)^-^^^^^  .  ^^^^)0.667 
a2  =  .  Cp^^gO'^Sa  .  k^^g^*^^^ 

^3  =  1^7,8°*''^^  •  Cp^^-^^^  •  ^40-667 

Subscripts  3  through  8  refer  to  refrigerant  key  locations  in  the  compressor  as 
per  figure  9.     If  separated  by  a  comma,  the  average  value  is  implied. 

The  heat  transfer  to/from  the  following  refrigerant  changes  the  refrigerant 
enthalpy  as  per  equation: 

q  =  mj.  '  M  (21) 

where       Ai  =  refrigerant  enthalpy  change 

Combining  equations  (18),  (19),  (20),  and  (21)  yields  the  following  expressions 
for  the  refrigerant  enthalpy  change  during  flow  between  respective  locations  of 
the  compressor: 

1     -  1     -  rn       .  k      0-667  .  p_      0.333  /t    _  t  \/fm  0-2  .  0.467v 

(22) 

i     -  i     -  rn        .V      0-667  .  p„      0.333  fr    -  r  ^/rm  0-2  .  ,,  0.466/ 
^6      ^7  -  ^^6,7      H,7  •  ^P6,7  ^^7      V^^'^r        '  ^6,7  > 

(23) 

-  '8  =  CQ6,7  (T7  -  l8)/(»r''''  •  > 

where        al,  a2,  a3  are  as  in  equation  (20) 

Derived  heat  transfer  relations  contain  heat  transfer  parameters  CQc  a>  ^Q4  5> 
CQ5  7,  and  CQ7  g.     These  parameters  are  functions  of  heat  transfer  surface  geo- 
metry mainly  and  have  to  be  found  empirically.     If  required  laboratory  test  data 
of  a  given  compressor  are  not  available,  typical  compressor  test  data  can  be 
used.     A  large  number  of  experimental  compressor  measurements  have  been  published. 
The  summary  of  these  data  and  a  list  of  references  can  be  found  in  [1] • 

Pressure  Drop  Relations 

Total  pressure  drop  AP^ot  experienced  by  a  flowing  fluid  results  from  pressure 
drops  due  to  friction,  momentum  change,  and  gravity,  i.e., 

^Ptot  =  ^Pfriction  +  ^^accel  +  ^Pgravity  (25) 

Pressure  drop  due  to  gravity  in  the  hermetic  compressor  may  be  disregarded 
based  on  an  order  of  magnitude  analysis.     On  the  same  grounds,  pressure  drop  of 
flowing  refrigerant  between  certain  compressor  locations  may  be  attributed  to 
either  viscous  effect  or  dynamic  effect. 

Pressure  drop  due  to  the  dynamic  effect,  i^'Paccel'        proportional  to  velocity 
head ,  i.e., 
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^^accel  -  P  •  V2  (26) 
which  can  be  written,  using  the  continuity  equation: 

^Paccel  =  CP  .  m//p  (27) 

where      .CP  =  pressure  drop  parameter 

nij.  =  refrigerant  mass  flow  rate 

V  =  refrigerant  velocity 

p  =  refrigerant  density 

The  relation  for  the  pressure  drop  due  to  viscous  effect,  Pfriction* 
derived  from  the  classical  Fanning  equation  for  pressure  drop  in  a  tube: 

^Pfriction  =  2f  .  p  .  V2  .  L/D  (28) 

where  f  =  friction  factor 

D  =  tube  diameter 
L  =  tube  length 

The  friction  factor,  f,  in  the  equation  above  is  approximately  proportional  to 
the  Reynolds  number  to  the  -0.2  power.  Considering  this  fact  and  applying  the 
equation  of  continuity,  equation  (28)  becomes: 

^Pfriction  =  CP  .  yO.2  .  ^^0.8/^  (29) 

Pressure  drop  of  the  refrigerant  in  the  compressor  can  is  modeled  by  evaluating 
individual  pressure  drops  between  refrigerant  key  locations  as  per  figure  9. 
Based  on  equations  (27)  and  (29),  and  attributing  pressure  drop  between  parti- 
cular compressor  locations  either  to  the  viscous  or  friction  effect,  the 
following  pressure  drop  relations  are  proposed: 

P^  -  P4  =  CP3  ^  •  m^.  /p^  4  for  compressor  can  inlet  (30) 

P4  -  P^  =  CP4  5  •        /P4  5  suction  manifold  and  (31) 

'  *  valve 

Pg  -  Py  =  CP^  -J  •  nij.  /p^  J  for  discharge  valve  and  (32) 

'  '  manifold 

Py  -  Pg  =  CPy  g  •  mj.-"-*^  •  My  8  dlscharge  line  (33) 

Pressure  drop  parameters  CP3^4,  CP4  5,  CP^  7,  and  CP7  3  have  to  be  found 
experimentally  for  a  given  compressor  or  can  be  calculated  using  test  data  for 
a  similar  compressor.    For  sources  of  such  data,  the  reader  may  refer  to  [1]. 

Equations  derived  in  this  section  make  it  possible  to  carry  out  an  energy 
balance  of  a  hermetic  compressor  solving  it  in  an  iterative  process  for  electri- 
cal energy  input,  for  heat  lost  to  the  ambient  air,  and  for  refrigerant  param- 
eters in  compressor  key  locations.     The  variety  of  designs  of  refrigerant  flow 
passages  in  the  compressor  caused  the  modeling  of  pressure  drop  and  heat  transfer 
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to  be  done  in  an  approximate  manner  in  this  general  compessor  model.     In  spite 
of  several  assumptions  that  were  made  to  simplify  the  modeling  process  and 
reduce  computing  time,  the  model  still  retains  sufficient  details  of  the  under- 
lying physical  principles  to  allow  designers  to  determine  which  specific  changes 
in  the  compressor  design  lead  to  increased  efficiency  of  the  compressor  and  the 
heat  pump  system. 

5.3     MODELING  OF  A  CONSTANT  FLOW  AREA  EXPANSION  DEVICE 

A  constant  flow  area  expansion  device,  as  used  in  heating/air  conditioning 
systems,  is  commonly  called  a  capillary  tube  or  a  refrigerant  flow  restrictor. 
Usually  it  is  a  small  bore  tube  of  length  as  short  as  one-half  inch  up  to  a  few 
feet,  connecting  the  outlet  of  the  condenser  (or  a  liquid  line)  to  the  inlet  of 
the  evaporator.     The  main  task,  of  the  constant  flow  area  expansion  device  is  to 
maintain  the  minimum  pressure  at  the  condenser  at  which  all  the  flowing  refri- 
gerant can  condense.    Many  researchers  have  investigated  flow  through  a  capil- 
lary tube  and  a  bibliography  on  the  subject  can  be  found  in  [4] .     The  objective 
of  this  research  was  to  explain  capillary  tube  phenomena  and  develop  a  mass 
flow  rate  prediction  method.     None  of  the  metods  found  in  the  literature  were 
computer  oriented  or  could  analyze  both  single-phase  and  two-phase  flow  in  the 
restrictor.     Thus  there  was  a  need  to  develop  a  more  comprehensive  capillary 
tube  routine  that  could  be  incorporated  into  a  heat  pump  simulation  program. 

5.3.1    Available  Capillary  Tube  Performance  Data 

Capillary  tube  is  a  traditionally  accepted  name  for  a  constant  flow  area 
expansion  device  used  in  heat  pump  systems.     This  name  is  adequate  and  mis- 
leading since  for  tube  diameters  in  the  neighborhood  of  1/16  of  an  inch, 
capillary  forces  are  negligible  for  the  task  the  capillary  tube  performs.  The 
pressure  drop  consists  practically  of: 

-  loss  due  to  sudden  contraction  at  the  entrance 

-  loss  due  to  flow  in  the  tube 

-  loss  due  to  sudden  enlargement  at  the  exit  to  the  evaporator 

The  flow  of  refrigerant  through  a  constant  bore  tube  is  more  complex  than  the 
geometric  simplicity  of  the  device  would  first  indicate.     The  pressure  and  tem- 
perature distribution  along  a  typical  capillary  tube  is  shown  in  figure  13. 
Bolstad's  description  of  the  flow  is  as  follows: 

At  the  entrance  to  the  tube,  section  0-1,  there  is  a  slight  pressure  drop 
which  was  usually  unreadable  on  the  gages .     From  point  1  to  point  2  the 
pressure  drop  is  linear.     In  the  portion  of  the  tube  0-1-2  the  refrigerant 
is  entirely  in  the  liquid  state  and  at  point  2  the  first  buble  of  vapor 
forms.     From  point  2  to  the  end  of  the  tube  the  pressure  drop  is  not 
linear,  the  pressure  drop  per  unit  length  increasing  as  the  end  of  the 
tube  is  approached.     For  this  portion  of  the  tube,  both  the  saturated 
liquid  and  saturated  vapor  phases  are  present,  the  percent  and  volume  of 
vapor  increasing  in  the  direction  of  flow  .... 
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With  a  saturation  temperature  scale  corresponding  to  the  pressure  scale 
superimposed  along  the  vertical  axis,  it  is  possible  to  plot  the  observed 
temperatures  in  a  more  meaningful  way  than  if  a  uniform  temperature  scale 
were  used.     The  temperature  is  constant  for  the  first  portion  of  the  tube 
0-1-2.     At  point  2,  the  presence  has  dropped  to  the  saturation  pressure 
corresponding  to  this  temperature.     Further  pressure  drop  beyond  point  2 
is  accompained  by  a  corresponding  drop  in  temperature,  the  temperature 
being  saturation  temperature  corresponding  to  the  pressure.     As  a  con- 
sequence, the  pressure  and  temperature  lines  coincide  from  point  2  to  the 
end  of  the  tube  [11] . 

The  point  of  the  tube  where  the  first  bubble  forms  is  called  a  bubble  point  or 
a  flash  point.    The  pressure  at  the  point  is  called  the  flash  pressure.  Bolstad 
also  presented  an  analytical  method  of  solution  for  adiabatic  flow  through  the 
capillary  based  on  the  Fanno  flow  theory. 

Mikol  [12]  performed  a  capillary  tube  investigation  from  which  his  conclusions 
can  be  summarized  as  follows: 

-  Fluid  flow  through  small  bore  tubes  conforms  to  continuum  flow  as 
established  for  large  bore  tubes  and  pipes. 

-  the  friction  factor  correlation  of  Moody  and  any  others  consistent 
with  Moody's  correlation  [13]  is  applicable  to  single-phase  flow  in 
small  bore  tubes. 

-  the  phenomena  of  metastability ,  persistence  of  the  liquid  state  to 
pressure  less  than  the  saturation  pressure  corresponding  to  its  local 
temperature,  has  been  found  to  occur. 

-  the  phenomenon  of  choked  flow  in  two-phase  flow  occurs  in  the  same 
way  and  for  the  same  reasons  as  in  the  case  of  gaseous  flow.  Sonic 
velocity  occurs  at  the  tube  exit. 

One  of  Mikol 's  findings,  existence  of  superheated  liquid  in  a  small  portion  of 
a  tube  was  not  observed  by  Bolstad  [11].     However,  it  was  reported  by  Cooper 
[14]  and  Rezk  [15] .     Investigators  have  found  that  delayed  evaporation  is 
affected  by  initial  disturbances  and  flow  agitation,  but  there  is  not  enough 
data  in  the  current  literature  to  assess  all  the  factors  promoting/eliminating 
this  phenomena  and  making  it  possible  to  consider  metastability  in  a  capillary 
tube  model  at  this  time. 


Refrigerant  mass  flow  rate  through  the  capillary  tube  always  increases  with  an 
increase  of  inlet  pressure.     It  also  increases  with  a  decrease  in  external 
outlet  pressure  down  to  a  critical  pressure  below  which  the  flow  does  not 
change.     Sensitivity  of  the  flow  to  external  outlet  pressure  decreases 
significantly  at  low  outlet  pressure  before  the  flow  is  choked. 

Pressure  and  temperature  distributions  along  a  capillary  tube  as  shown  in  figure 
13  occur  at  design  operating  conditions  of  a  long  (a  few  feet)  capillary.  Part 
of  the  capillary  is  filled  with  flowing  liquid,  while  two-phase  flow  exists  in 
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DISTANCE  FROM  CAPILLARY  INLET  (Ft) 


Figure  13.    Pressure  and  temperature  distribution  along  typical 
capillary  tube  [11] 
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the  other  part.     However,  there  are  also  other  possible  modes  of  operation, 
i.e.,  with  only  two-phase  flow  in  the  capillary  (the  case  of  Incomplete  conden- 
sation in  the  condenser)  or  with  only  liquid  flow  (the  case  of  a  short  restric- 
tor).     All  these  cases  are  observed  in  practice  and  have  to  be  simulated  by  a 
general  model  of  the  constant  flow  area  expansion  device. 

Based  on  the  experimental  evidence  and  the  theory  of  large  tube  fluid  mechanics, 
the  following  assumption  were  set  up  for  model  formulation: 

1.  The  capillary  tube  is  a  straight,  horizontal,  constant  inner  diameter 
tube 

2.  Flow  in  the  capillary  is  one-dimensional  and  homogeneous 

3.  Flow  in  the  capillary  is  adiabatic 

4.  Flow  resistance  in  the  capillary  tube  can  be  subdivided  into 

a.  resistance  due  to  the  entrance  effect 

b.  resistance  due  to  flow  in  a  tube  which  consists,  in  the  general 
case,  of  single-phase  liquid  flow  resistance  from  the  entrance 
to  the  flash  point,  and  two-phase  mixture  flow  resistance  in  the 
rest  of  the  tube.     Existence  of  delayed  evaporation  phenomena 

is  neglected.  Resistance  due  to  the  exit  effect  is  neglected  as 
meaningless  for  a  choked  flow  and  insignificant  for  a  non-choked 
flow  [11] . 

5.  Choked  flow  phenomenon  for  two-phase  is  governed  by  the  same  laws  as 
for  single-phase  flow,  and  can  be  modeled  accordingly. 

5.3.2     Critical  Flow  in  a  Capillary  Tube 

Refrigerant  mass  flow  rate  through  a  given  capillary  tube  will  depend  on  the 
inlet  refrigerant  state  and  on  pressure  that  will  be  established  at  the  tube 
outlet.     This  pressure  may  equal  the  evaporator  pressure  or  may  have  a  higher 
value  if  the  flow  is  choked.     Since  pressure  at  the  capillary  tube  exit  is  one 
of  the  parameters  affecting  the  flow,  it  has  to  be  known  for  proper  refrigerant 
mass  flow  rate  prediction. 

The  assumption  that  the  flow  in  the  capillary  is  one-dimensional  and  homogeneous 
enables  the  two-phase  flow  in  the  capillary  to  be  treated  as  single-phase  flow 
with  uniform  properties  at  any  cross-section  of  the  flow  and  allows  use  of  the 
single-phase,  one-dimensional  form  of  the  governing  equations  as  presented  in 
figure  14,  15,  and  16. 

Adiabatic  flow  through  a  capillary  tube  is  a  classical  example  of  so  called 
Fanno  flow;  adiabatic  flow  with  friction  in  a  constant  area  duct.     The  energy 
equation  for  such  flow  has  the  form: 

di  +  V  •  dV  =  0        or        1^=1+  V^/2  =  const  (34) 
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which,  when  combined  with  the  equation  of  continuity  yields: 


io  =  i  +  —  •  v^/2  =  i  +        •  v^/2  (35) 

where  A  =  tube  cross-sectional  area 

m 

G  = — E-  ,  refrigerant  mass  flux 
A 

i  =  refrigerant  enthalpy 
io  =  refrigerant  stagnation  enthalpy 
THj.  =  refrigerant  mass  flow  rate 

V  =  flow  velocity 

V  =  refrigerant  specific  volume 

Graphical  representation  of  equation  (35)  on  the  enthalpy-entropy  diagram,  as 
shown  in  figure  17,  is  called  a  Fanno  line.     Fanno  flow  as  an  irreversible 
adiabatic  process,  can  sequentially  exist  only  in  the  direction  of  increasing 
entropy.     The  upper  branch  of  the  Fanno  line  corresponds  to  subsonic,  acceler- 
ating flow  while  the  lower  branch  is  the  supersonic,  decelerating  flow.  Both 
flows  tend  towards  point  C  where  the  sonic  velocity  is  reached. 

The  Fanno  line  implies  that  fluid  cannot  reach  the  velocity  of  sound  inside  a 
constant  area  duct,  because,  if  this  happened,  further  flow  in  the  duct  would 
have  to  be  associated  with  a  decrease  of  entropy,  as  per  figure  17,  and  that 
would  be  in  violation  of  the  Second  Law  of  Thermodynamics.     Thus,  the  only 
location  of  a  tube  at  which  sonic  velocity  can  be  reached  is  the  tube  exit,  and 
choking  will  never  occur  inside  the  capillary  regardless  of  external  outlet 
pressure.     Evaluation  of  the  critical  pressure  for  two-phase  Fanno  flow  is 
discussed  in  appendix  E. 

It  is  intuitively  felt,  based  on  compressible  gas  theory,  that  if  the  critical 
pressure  is  greater  than  the  pressure  in  an  evaporator,  the  critical  pressure 
will  be  established  at  the  capillary  tube  outlet.     However,  it  is  not  clear 
what  pressure  would  actually  exist  at  the  end  of  a  capillary  which  is  fed  with 
a  subcooled  liquid,  if  the  calculation  choking  pressure  is  higher  than  the  sat- 
urated liquid  pressure  resulting  from  the  liquid  temperature^^/ .     Before  this  is 
explained,  some  information  on  sonic  velocity  phenomena  is  reviewed. 

The  velocity  of  sound  is  calculated  by  the  formula: 
0.5 

(36) 

s 


3/ 

It  is  assumed  here,  for  discussion  purposes  only,  that  choking  pressure  can 
be  found  independently  from  refrigerant  mass  flow  rate  and  then  mass  flow 
rate  can  be  calculated.     In  fact,  mass  flow  rate  and  choking  pressure  (and 
other  factors  of  operation)  are  interrelated  and  have  to  be  found 
simultaneously  in  an  iterative  process. 
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Figure  14.    Mass  balance  for  an  element  of  fluid  in  one-dimensional 
steady  flow  in  a  constant  area  duct 
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Figure  15.    Momentum  balance  for  an  element  of  fluid  in  one-dimensional 
steady  flow  in  a  horizontal,  constant  area  duct 
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Figure  16.     Energy  balance  for  an  element  of  fluid  in  one-dimensional 
adiabatic,  steady  flow  in  a  horizontal,  constant  area  duct 
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Figure  17 .     Fanno  line 
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For  a  liquid,  in  which  the  density  is  hardly  affected  by  pressure,  the  velocity 
of  sound  assumes  a  very  high  value,  typically  in  the  range  of  5000  ft/s.  Liquid 
flow  in  capillary  tubes  is  never  accelerated  to  such  high  velocities.  However, 
the  velocity  of  sound  for  a  two-phase  state  containing  a  small  amount  of  vapor 
is  drastically  smaller,  e.g.,  sonic  velocity  of  two-phase  flow  of  quality  less 
than  0.0001  and  temperature  of  100°F  is  on  the  order  of  1  ft/s.     Sonic  velocity 
increases  with  temperature  and  quality  and  assumes  a  value  of  about  500  ft/s 
close  to  the  critical  state.     Thus,  the  velocity  of  sound  in  a  liquid  is  much 
different  than  in  the  two-phase  state  and  a  sharp  discontinuity  in  sonic 
velocity  takes  place  at  the  boundary  of  single  and  two-phase  flow  [16] . 

In  the  case  of  a  capillary  tube  working  at  design  conditions,  there  is  a  single- 
phase  liquid  flow  in  the  beginning  part  of  the  tube  and  single  component  two- 
phase  mixture  flow  in  the  latter  tube  portion.     Though  there  is  basically  one 
substance  flowing,  the  relation  between  pressure  and  density  along  the  path  of 
flow  will  be  different  for  each  part  of  the  flow  in  the  tube.     Theoretically,  a 
separate  Fanno  line  could  be  drawn  for  each  flow,  with  an  intersection  point 
corresponding  to  a  flash  point.     Also,  two  sonic  velocities  could  be  considered, 
but  a  limitation  of  sonic  velocity  in  the  liquid  is  disregarded  as  it  is  too 
high  to  be  reached  in  a  capillary  tube  under  any  operational  conditions. 

Consider  several  capillary  tubes  of  the  same  bore  but  different  lengths.  All 
the  capillaries  receive  subcooled  liquid  at  the  same  condition  and  are  working 
at  the  same  pressure  in  the  evaporator.     The  various  possible  flows  are  shown 
in  figure  18.     If  a  capillary  tube  is  sufficiently  long  (a),  it  will  be  able  to 
realize  a  pressure  drop  down  to  the  evaporator  pressure  and  choking  will  not 
occur.     With  a  shorter  capillary  tube  such  as  (b)  and  (c),  the  mass  flow  rate 
is  greater,  and  the  critical  pressure  which  is  established  at  the  tube  exit, 
exceeds  the  pressure  in  the  evaporator.     Tube  (d)  has  the  length  at  which  the 
critical  pressure  equals  the  flash  pressure.     The  critical  pressure  in  the  case 
of  tube  (e)  is  greater  than  the  flash  pressure  and  the  latter  is  established  at 
the  tube  exit  due  to  the  flash  that  will  occur  from  the  sudden  expansion. 
Thus,  any  tube  shorter  than  tube  (d)  will  have  flash  pressure  at  the  exit.  To 
prove  this  last  case,  which  is  least  evident,  other  flow  situations  than  that 
proposed  must  be  considered.     First,  eliminated  the  possibility  of  pressure  at 
the  tube  exit  being  greater  than  the  flash  pressure  because  otherwise  it  would 
imply  the  liquid  is  choked  at  other  than  its  critical  pressure,  which  is 
physically  impossible.     Another  alternative  is  to  assume  that  the  liquid  reaches 
saturation  pressure  within  the  tube  and  changes  into  a  two-phase  mixture.  With 
the  two-phase  mixture  critical  pressure  being  higher  than  saturation  pressure, 
it  would  mean  a  change  from  single-phase  subsonic  flow  into  two-phase  super- 
sonic flow.     Such  flow  is  also  impossible  since  it  would  require  the  flow  to 
experience  a  decrease  in  quality  in  its  progress  and  would  be  converted  back  to 
liquid  flow  immediately.     The  liquid  flow  also  could  not  proceed  any  further 
because  vapor  bubbles  would  form.     Thus  the  only  possible  situation  for  flow 
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CAPILLARY  TUBE  LENGTH 

*  Actual  flash  pressure  is  not  constant.  It  Increases  slightly  with 
capillary  tube  length. 

Figure  18.    Pressure  distribution  along  different  length  capillary 
tubes  at  the  same  operating  conditions 
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with  the  critical  pressure  greater  than  the  liquid  saturation  pressure  is  to  be 
choked  at  the  tube  exit  where  the  flow  is  at  the  liquid  saturation  pressure. 

5.3.3    Model  Formulation 

Refrigerant  flow  from  the  liquid  line  into  the  capillary  tube  experiences  a 
pressure  drop  due  to  sudden  contraction.     This  pressure  drop  consists  of  an 
acceleration  loss  and  entrance  friction  loss  and  is  usually  expressed  by  a 
decrease  in  the  Bernoulli  head  and  a  contraction  coefficient  referred  to  the 
kinetic  energy  of  the  flow  in  the  section  of  smaller  flow  area: 


(37) 


Subscripts  in  eq.  (37)  refer  to  sections  shown  in  figure  13.     Combining  with 
the  equation  of  continuity,  the  above  can  be  rearranged  to: 

^1^ 

Pq  -        =  AP  =  (1  +  K)Po^i  '-j-  (38) 


The  contraction  coefficient,  K,  given  in  the  literature  is  stricly  empirical 
and  is  represented  as  a  function  of  the  contraction  area  ratio.     Several  sources 
are  in  disagreement  about  its  value.     The  value  of  K  =  0.15,  used  here  is  from 
a  derivation  based  on  Kay's  general  formula  [17].     It  was  calculated  for  a 
normal  range  of  contraction  area  ratios  for  capillary  tubes  with  slightly 
beveled  entrances. 

The  equation  of  motion  for  steady  flow  in  a  constant  cross-section  area  pipe 
has  the  following  form: 

p'A'dV+A'dP+T«S»dL=0  (39) 

The  skin  friction  coefficient,  x,  can  be  expressed  in  terms  of  the  friction 
coefficient,  f,  and  the  velocity  head: 

T  =  1  f  .  p  .  v2  (40) 
2 

The  flow  velocity  term,  V,  can  be  eliminated  by  means  of  the  equation  of 
continuity: 

d(V  •  p)  =  0  (41) 
Substituting  and  rearranging;  the  equation  of  motion  assumes  the  following  form: 

2      ^i+1  ^i+1 

^      •  /    p  •  dP  +  1  •  /     f  •  dL  +  In  -^-i —  =  0  (42) 
mr  D  p^+i 

Pi  Li 
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As  discussed  before,  flow  in  a  capillary  tube,  in  the  general  case,  can  be 
subdivided  into  two  parts  separated  by  a  flash  point:     the  liquid  flow  part  and 
the  two-phase  mixture  part.     The  same  equation  will  be  applicable  for  both 
flows  through  in  the  case  of  liquid  flow,  it  can  be  simplified  on  grounds  of 
incompressibility .     In  fact,  it  reduces  to  the  Fanning  pressure  drop  formula  in 
the  following  form: 


AP  =^ 


g2  .  L 


(43) 


where  the  friction  factor,  f,  can  be  evaluated  by  the  empirical  formula: 

f  =  0.046  •  Re"0-2  (44) 

For  the  two-phase  mixture  flow,  equation  (42)  has  to  be  solved  in  its  full 
form.     This  was  done  by  Whitesel,   [18]  and  [19],  but  with  significant  over- 
simplified approximations  for  refrigerant  properties.     In  solving  equation  (42), 
difficulty  arises  with  evaluating  the  first  term  because  it  depends  directly  on 
the  pressure-density  relation  along  the  path  of  flow.     However,  the  relation 
can  be  obtained  by  considering  the  adiabatic  flow  case.     The  specific  volume  at 
a  given  pressure  can  be  expressed  in  terms  of  the  property  values  for  saturated 
liquid  and  vapor  and  in  terms  of  quality: 

V  =  vl  +  x(Vv  -  vl)  (45) 

where         v  =  specific  volume 
X  =  quality 

Subscripts  L  and  V  are  for  liquid  and  vapor,  respectively 

The  quality  of  the  flow  in  the  Fanno  path  can  be  found  as  explained  in  appendix 
E.     Thus  integration  of  refrigerant  density  over  a  given  pressure  interval  can 
be  done  numerically.     Still  another  problem  is  faced  in  evaluating  the  second 
term  of  equation  (42),  which  includes  a  two-phase  friction  factor  as  a  function 
of  tube  length.     Erth  [20]  made  an  effort  to  correlate  two-phase  average  fric- 
tion factor  in  a  capillary  tube.     His  regression  analysis,  based  on  four  sets 
of  data  from  four  different  experiments,  yielded  the  following  correlation  for 
a  two-phase  flow  mean  friction  factor,  f^,  as  a  function  of  the  inlet  conditions 
only: 

=-^^exp[(l  -  xi        )/2.4]  (46) 
Re^*^ 

where        x^  =  quality  of  refrigerant  entering  capillary  tube 

Re  =   ^  *  P   ,  two-phase  Reynolds  number  (47) 

PL     ^i  (^^V  ~  ^l) 

Using  this  mean  friction  factor,  f^j,  equation  (42)  may  be  written  for  the  two- 
phase  portion  of  the  flow  in  the  following  form: 
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0.5 


mj.  =  A 


I  f„  •  (L3  -  4)  +  m  (P2/P3) 


(48) 


where  '  numbers  used  as  subscripts  denote  location  consistent  with  figure  13. 

The  above  equations  have  to  be  solved  in  a  highly  iterative  process  since 
choking  pressure,  friction  factor,  fraction  of  capillary  tube  length  with  liquid 
and  two-phase  flow,  and  the  velocity  head  used  to  correct  enthlapy  are  functions 
of  refrigerant  mass  flow  rate  which  has  yet  to  be  found. 

5.4     MODELING  OF  AN  EVAPORATOR  AND  A  CONDENSER 

5.4.1    Modeling  Methodology 

There  are  two  heat  exchangers  in  a  heat  pump:     an  indoor  coil  and  an  outdoor 
coil.     Both  coils  are  made  in  a  similar  way  and  both  serve  as  an  evaporator  or 
condenser  depending  on  the  heat  pump  opertional  mode.     Schematic  of  a  typical 
heat  pump  heat  exchanger  is  shown  in  figure  19.     It  consists  of  a  set  of  finned 
tubes  connected  in  a  specifically  designed  circuit  configuration..     An  example 
of  the  coil  circuitry  is  illustrated  in  figure  20.     The  refrigerant  is  pumped 
through  the  tubing  while  air  is  forced  to  flow  over  the  outside  of  the  coil. 

Heat  exchanger  models  in  the  MIT  and  ORNL  heat  pump  simulation  programs,  [1] 
and  [2],  are  very  similar.     Coil  performance  calculations  in  these  models  are 
based  on  the  ef f ectiveness-NTU  method  for  a  cross-flow  heat  exchanger  with  both 
fluids  umixed.     These  models  can  rigorously  simulate  only  certain  regular  tube 
arrangements  and  are  Insenltive  to  different  coil  circuitry  patterns.     The  heat 
exchanger  models  in  [3]  are  based  on  the  so  called  tube-by-tube  simulation 
method  which  depends  on  imaginary  isolation  of  one  tube  with  appropriate  fin 
surfaces  from  the  coil  assembly  and  calculating  the  performance  of  this  tube 
independent lyA/ .     The  heat  transfer  to/from  a  tube  is  calculated  with  the  aid 
of  heat  exchanger  cross-flow  theory.     This  is  done  for  each  tube  in  a  proper 
sequence  and  the  summation  results  in  total  coil  capacity.     An  even  more  detailed 
method  of  heat  exchanger  simulation  would  be  to  apply  the  finite-difference 
technique.     Such  a  model  would  subdivide  an  analyzed  tube  of  the  coil  into 
arbitrary  small  elements  and  would  solve  heat  transfer  equations  for  each  ele- 
ment marching  in  space  from  the  tube  inlet  to  the  tube  outlet.     Since  for  this 
study  the  cost  of  computing  time  required  for  such  simulation  would  be 
prohibitive,  the  tube-by-tube  approach  was  adopted. 

In  the  tube-by-tube  simulation  method,  the  performance  of  one  tube  is  calculated 
at  one  time  and  the  general  finned  coil  problem  is  reduced  to  a  single  tube  pro- 
blem.    In  order  to  conduct  heat  transfer  calculations,  fin  surfaces  associated 


ZJ  At  the  stage  of  completion  of  this  report  a  new  tube-by-tube  computer  model 
of  a  condenser  coil  was  published  by  ORNL  [56] . 


50 


Figure  19.     Schematic  of  a  heat  pump  heat  exchanger 
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Figure  20 .     Example  of  coil  circuitry 
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with  the  tubes  must  be  defined.     Following  Carrier  and  Anderson  [22],  It  was 
assumed  that  fin  area  served  by  each  tube  Is  equivalent  In  performance  to  a 
circular-plate  fin  of  equal  area.     Thus  a  single  tube  is  considered  with  a 
circular  fin  of  diameter,  D^,  as  shown  in  figure  21. 

5.4.2    Heat  Transfer  Rate  for  a  Tube  in  a  Cross-Flow  Arrangement 

Usually  a  heat  pump  coll  employs  some  form  of  cross-flow  arrangement.     If  a 
separate  tube  Is  considered,  the  problem  is  one  of  pure  cross-flow.  Fortunately, 
this  kind  of  arrangement  has  received  much  attention  in  theoretical  investiga- 
tions.   According  to  the  Peclet  heat  transfer  equation: 


Q  =  U  •  Ah  •  AT 


(49) 


where       Aj^  =  heat  transfer  surface  area 

U  =  overall  coefficient  of  heat  transfer 
AT  =  temperature  difference 

In  case  of  a  pure  cross-flow  arrangement,  one  of  two  mean  temperature  differences 
applies  for  AT,  [23]: 


t2  -  ti 


To  -  ti 

In^^  1- 


Ti  -  t- 


when  temperature  of  one  of  the 
fluids  does  not  change  (tem- 
perature of  that  fluid  is 
denoted  by  T) 


(50) 


m 


In 


t2  -  ti 


Tl  -  T2 
^2  -  ^1 


T,   -  T,  T2  _  t, 
 ^+  In    ^  ^ 


t2  -  ti 


Ti  -  ti 


when  temperature  of  both 
fluids  changes 


(51) 


where         T  =  temperature  of  one  fluid 

t  =  temperature  of  another  fluid 
ATjQ  =  mean  temperature  difference 

subscripts  1  and  2  refer  to  tube  inlet  (upstream)  and  outlet 
(downstream)  conditions 

The  enthalpy  change  of  the  fluid  will  be  according  to  the  equation: 
Q  =  m  •  (12  -  ii) 

or 

Q  =  m  •  Cp  •  (T2  -  Ti)  (52) 
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Figure  21.  Approximation  method  for  treating  a  rectangular-plate 
fin  of  uniform  thickness  in  terms  of  a  flat  circular- 
plate  fin  of  equal  area 


54 


where        Cp  =  specific  heat  of  fluid 
i  =  enthalpy 
m  =  mass  flow  rate 


Looking  at  any  vapor  compression  cycle  P-h  or  T-s  diagram,  it  can  be  realized 
that  the  single-phase  and  two-phase  refrigerant  flow  usually  exists  in  a  given 
heat  exchanger.     Also,  both  flow  patterns  can  actually  exist  in  one  tube.  That 
means  the  refrigerant  flowing  in  the  tube  may  maintain  approximately  constant 
temperature  (two-phase)  or  its  temperature  may  vary  (single-phase).     Not  only 
the  mean  temperature  difference  between  fluids  but  also  the  refrigerant  pressure 
drop  is  affected  by  the  flow  pattern.     Equations  presented  below,  derived  from 
equations  (49),  (50),  (51),  and  (52),  allow  for  detailed  consideration  of  these 
problems.     The  heat  transfer  rate  for  each  mentioned  flow  condition  can  be 
calculated  as  follows: 

-  single-phase  flow  only,  refrigerant  is  superheated  or  subcooled  at  both 
inlet  and  outlet 

Q  =  Cpr  •  mr(T^^i  -  Ta,i)(l  -  exp(  -  \        (1  -  exp(-  ^  '  ))))  (53) 

Pr  *     r  Pa  *  a 

-  two-phase  flow  only,  refrigerant  is  in  two-phase  at  both  inlet  and 
outlet 

Q  =  Cpa  .        (1  -  exp(  '  ))(Tr,g  "  Ta,i)  (54) 

"a  a 

-  superheated  vapor  at  tube  inlet,  two-phase  at  the  tube  outlet 


U  . 

Q  =  m^-dr  i  -  ir  V)  +  Cpa  •  ma(l  -  Zy)  (1  -  exp(-  ^))(Tr,g  -        i)  (55) 

where  =  fraction  of  the  tube  length  in  the  superehated  region  which 

can  be  calculated  by  the  equation: 

-''^ '  • '''''  -  cp'(<  [ 

z   =  L^ii  (56) 

V  U   .  Aq 

CPa  •  mad  "  exp(l  -  —  — )) 

^Pa  •  "'a 

-  two-phase  at  tube  inlet,  subcooled  liquid  at  the  tube  outlet 

Cp^  .  m,(l-Z^„)  U  .  A- 

Q  =  Cp,  .  .,(T,,g  -  T,,i)(l  -  axp(-    \  (1  -  exp(-^^_^)))) 

r  *^a  a 

+  m^d^^i  -  i^^L)  (57) 
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where      Zj-p  =  fraction  of  the  tube  length  in  the  two-phase  region  which  can 
be  calculated  by  the  equation: 


 '°r^^r,i  "  ^r,L^  

Ztp  =  U  •  A-  (58) 

Cpa  •  mad  -  exp(-  ^p    .  I  ^^^'^r.i-  Ta,i) 

In  equations  (53)  through  (58)  the  following  nomenclature  was  used: 


Aq  =  total  exterior  surface  area  associated  with  the  tube  wetted  by  air 
Cpa  =  air  specific  heat  at  constant  pressure 
Cpj.  =  refrigerant  specific  heat  at  constant  pressure 
^r  i  ~  refrigerant  enthalpy  at  the  tube  inlet 
ij.  L  =  enthalpy  of  refrigerant  saturated  liquid 
ij.^y  =  enthalpy  of  refrigerant  saturated  vapor 

=  air  mass  flow  rate  associated  with  the  tube 
rnj.  =  refrigerant  mass  flow  rate  in  the  tube 
£  =  air  temperature  upstream  of  the  tube 
Tj-^g  =  refrigerant  saturation  temperature 
Tj.  i  =  refrigerant  temperature  at  the  tube  inlet 
U  =  overall  tube  heat  transfer  coefficient 


A  desirable  feature  in  the  above  equations  is  that  they  do  not  involve  any  of 
the  outlet  temperatures  or  any  of  the  outlet  fluid  parameters.    Hence,  the  heat 
transfer  rate  can  be  calculated  based  on  total  heat  transfer  area,  overall  tube 
heat  transfer  coefficient,  inlet  fluid  parameters,  and  their  mass  flow  rates. 

5.4.3    Refrigerant  and  Air  Mass  Flow  Rates  Associated  with  a  Tube 


Refrigerant  Mass  Flow  Rate 


During  flow  through  a  heat  pump  heat  exchanger,  the  refrigerant  undergoes  a 
change  of  phase  in  the  course  of  evaporation  or  condensation.     The  change  of 
phase  is  associated  with  a  dramatic  change  of  density  which  affects  the  velo- 
city and  pressure  drop  of  the  working  fluid.     In  order  to  prevent  a  high  pres- 
sure drop,  tubes  in  some  heat  pump  coils  are  connected  to  form  branched  circuits. 
An  example  of  such  coil  circuitry  has  been  shown  in  figure  20. 

Refrigerant  flow  direction  marked  in  figure  20  is  for  the  coil  working  as  a 
condenser.     Superheated  vapor  enters  the  vapor  side  manifold  and  is  distributed 
into  10  tubes.     On  its  flow  path,  refrigerant  merges  several  times  and  finally 
merges  in  the  liquid  manifold  to  enter  the  liquid  line.     For  the  coil  operating 
as  an  evaporator,  the  direction  of  flow  is  opposite  to  that  marked  in  the 
figure.    Low  quality  refrigerant  enters  the  coll  and  the  flow  is  subdivided 
into  three  circuits.     On  its  way  through  the  coil  the  refrigerant  evaporates 
and  splits  several  times  on  its  way  towards  the  exist  where  it  is  finally  col- 
lected into  one  bigger  diameter  vapor  line.     The  mass  flow  rates  through  the 
particular  circuits  of  the  coil  are  self-adjusting  so  the  pressure  at  merging 
(splitting)  tubes  is  the  same. 
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To  perform  a  simulation  of  a  coil  by  the  tube-by-tube  method,  refrigerant  mass 
flow  rate  for  each  tube  has  to  be  known.     Since  total  refrigerant  mass  flow 
rate  supplied  to  the  coil  is  known,   the  problem  reduces  to  the  determination  of 
refrigerant  distribution.     In  [3],  refrigerant  distribution  is  found  by  consid- 
ering the  coil  working  as  an  evaporator  and  assuming  that  refrigerant  splits 
equally  among  tubes  at  any  splitting  point  of  tubing  circuitry.     Thus  the  frac- 
tion of  the  total  flow  is  assigned  to  each  tube  based  on  circuit  geometry  and 
is  used  in  both  the  heating  and  cooling  mode  simulation. 

The  method  for  determination  of  refrigerant  flow  distribution  used  in  [3]  was 
checked  by  examining  calculated  refrigerant  pressures  at  the  ends  of  all  cir- 
cuit branches.     These  pressures  should  be  equal  if  the  assumed  distribution  is 
right;  however,  the  pressures  were  found  to  be  different  by  as  much  as  3  psi. 
Of  course,  refrigerant  flow  distribution  could  be  corrected  in  an  iteration 
process.     This  would  require  repeating  all  coil  simulation  calculations  several 
times,  hence  making  coil  simulation  prohibitively  expensive. 

Another  approach  was  tried  to  determine  refrigerant  flow  distribution.  Since 
most  of  the  coil  total  pressure  drop  may  be  expected  to  result  from  superheat 
vapor  and  two-phase  flow,  it  was  assumed  that  the  refrigerant  flow  is  uniformly 
distributed  among  tubes  connected  to  the  vapor  side  manifold  and  that  mass  flow 
rates  in  other  tubes  may  be  found  by  following  the  refrigerant  path  with  direc- 
tion of  flow  as  marked  in  figure  20.     The  resulting  distribution  was  checked  as 
before.     The  maximum  pressure  discrepancy  found  was  equal  to  0.3  psi  which 
represented  less  than  a  0.2°F  variation  in  the  saturation  temperature  between 
merging  tubes.     This  was  considered  satisfactory  and  the  method  was  adopted  in 
this  study  for  determining  refrigerant  mass  flow  rate  distribution  in  a  coil. 


Air  Mass  Flow  Rate 


Air  mass  flow  rate  was  assumed  to  be  distributed  uniformly  over  the  whole 
coil  face  regardless  of  the  coil  and  fan  respective  locations,  so  each  coil  was 
associated  with  the  same  air  mass  flow.     It  was  also  assumed  that  there  were 
enough  turbulence  in  the  air  stream  passing  through  the  coil  to  provide 
effective  mixing  so  that  air  of  uniform  properties  enters  each  tube  bank. 

5.4.4    Overall  Heat  Transfer  Coefficient  for  a  Dry  Finned  Tube 

Dry  fined  tube  analysis  is  applicable  to  a  condenser  and  also  to  an  evaporator 
if  no  dehumidif cation  takes  place.     The  overall  heat  transfer  coefficient,  U, 
for  a  dry  finned  tube  can  be  derived  by  summing  the  individual  resistances 
between  the  refrigerant    and  the  air,  [23]: 


U  = 


Ap.ihi 


Ap,ihd,i 


o  P 
Ap , m^p 


^(1  - 
c,o^ 


Af 


(1  -  <{))) 


(59) 


where       Af  =  fin  surface  area 

Aq  =  total  exterior  surface  area  wetted  by  air 
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Ap^£  =  pipe  inside  surface  area 
^p,m  ~  PiP^  mean  surface  area 

h(,  Q  =  convection  heat  transfer  coefficient  at  exterior  surface 
i  ~  heat  transfer  coefficient  for  inside  tube  deposit 
Aj^  =  inside  tube  convection  heat  transfer  coefficient 
kp  =  thermal  conductivity  of  pipe  material 
Xp  =  thickness  of  pipe  wall 

T        —  T 
(j)  =  -J^  _a_  ^  fin  efficiency 


Tf.b 


-  T 


a 


=  air  temperature 
Tf      =  fin  base  temperature 
^f,m  ~  niean  fin  temperature 

The  third  term  of  equation  (59)  can  be  evaluated  if  the  heat  exchanger  material 
and  geometry  are  known.  The  heat  transfer  coefficient  for  inside  tube  deposit, 
h^  i>  assumed  to  be  5000  Btu/(h  'ft     •  F)  [23]  which  allows  for  evaluation 

of 'the  second  term.     Terms  1  and  4,  which  refer  to  the  inside  and  outside  con- 
vection resistance  respectively,  required  considerably  more  analysis  to 
establish  the  proper  algorithm  for  determining  the  heat  transfer  coefficient. 

5.4.5    Forced  Convection  Heat  Transfer  Inside  a  Tube 

Analyzing  the  problem  for  both  an  evaporator  and  a  condenser,  the  following 
modes  of  forced  convection  are  encountered: 

-  single-phase  forced  convection 

-  two-phase  forced  convection  with  condensation 

-  two-phase  forced  convection  with  evaporation 

The  physics  of  these  phenomena  are  very  much  different  and  they  have  to  be 
considered  separately. 

Single-Phase  Forced  Convection 

Single-phase  forced  convection  takes  place  in  a  condenser,  at  the  entrance 
section  where  the  superheated  vapor  is  being  cooled,  and  a  the  exit  section 
where  a  subcooled  liquid  is  being  cooled.     It  is  also  applicable  in  the  eva- 
porator, as  the  superheated  vapor  passes  through  the  exit  tubes.     The  non- 
dimensional  heat  transfer  parameter  describing  this  phenomena,  Nusselt  number, 
is  related  to  the  non-dimensional  Reynolds  and  Prandtl  numbers  in  the  following 
form: 

Nu  =  0.023  •  ReO-8  .  PrO.333  (60) 

where  Nu  =  ^  *  ^,  Nusselt  number 

k 

Pr  =iL_!_C£,  Prandtl  number 
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Re  s  G  •        Reynolds  number 

y 

Cp  =  specific  heat  at  constant  pressure 

D  =  Inside  diameter  of  a  tube 

G  =  refrigerant  mass  flux 

h  =  convection  heat  transfer  coefficient 

k  =  refrigerant  thermal  conductivity 

y  =  refrigerant  absolute  viscosity 

Two-Phase  Forced  Convection  with  Condensation 

The  predominant  flow  pattern  during  condensation  in  a  heat  pump  condenser  is 
annular  flow  with  liquid  refrigerant  flowing  on  the  pipe  wall  and  vapor  refri- 
gerant flowing  in  the  core.     Traviss,  Baron,  and  Rohsennow  [24]  performed  experi- 
mental and  analytical  studies  of  vapor  condensation  in  annular  flow  with 
refrigerants  R-12  and  R-22.     They  assumed  von  Karman  universal  velocity  distri- 
bution in  the  condensate  flow  (like  on  a  flat  plate),  calculated  pressure  drops 
using  the  Lockhart-Martinelli  method  [25],  and  applied  momentum  and  heat  transfer 
analogy  for  calculation  of  the  condensation  heat  transfer  coefficient.  The 
correlation  proposed  by  them  has  the  following  form: 

Re.^*^  •  Pr.  •  Fl^ 

Nu  = —  ^   (61) 

F2 

where        jj^  =  h  «  D 

h  =  condensation  heat  transfer  coefficient 
D  =  tube  Inside  diameter 
kx,  =  thermal  conductivity  of  liquid  refrigerant 

R      =  G(l  -  x)D 

G  =  refrigerant  mass  flux 
x  =  quality 
HL  -  liquid  refrigerant  absolute  viscosity 

L 

3  =  1  for  Fl  <  1,  3  =1.15  for  Fl  >  1 
Fl  and  F2  in  equation  (61)  are  dimensionless  parameters  expressed  as  follows: 

Fl  =  0.15  (X^^-1  +  2.85  X,.^^-^24)  (62) 

F2  =  0.707  •  Pr^  •  ^^j^^'^  ^^L  <  5° 
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F2  =  5  •  Pr^  +  5  •  ln(l  +  PrL(0. 09636  •  for  50  <  Re^  <  1125 


P2  =  5 


PrL  +  5 


ln(l  +  Pr^)  +  2.5  •  ln(0. 00313  •  Re^ 


0.812 


)  (63) 
for  ReL  <  1125 


Parameter  ^x:t>  formulated  by  Lockhart-Martinelli  [25]  with  the  assumption  of  no 
radial  pressure  gradient  and  a  smooth  pipe,  has  the  following  form: 


"■tt 


1  -  X 


0.9 


0.5 


0.1 


(64) 


Parameter  X^t  Is  inversely  proportional  to  flow  quality  and  refers  to  turbulent 
vapor  and  turbulent  liquid  flow.     Physically  it  is  equal  to  the  square  roof  of 
the  ratio  of  pressure  drop  of  the  liquid  phase  to  the  frictional  pressure  drop 
of  the  vapor  phase  if  these  phases  were  flowing  independently  and  alone  in  the 
tube ,  i.e.: 


tt 


dP 
dL 


dP 

dL  /V 


0.5 


(65) 


Equation  (61)  is  applicable  where  conditions  for  annular  condensation  in  a 
tube  exist.     Such  conditons  may  be  assumed  to  exist  for  flow  qualities  ranging 
from  0.1  to  0.9.    At  qualities  larger  than  0.9,  the  whole  tube  inner  surface 
is  not  covered  by  a  liquid  film  and  part  of  the  heat  transfer  is  just  that  of 
single-phase  convection.    At  qualities  less  than  0.1,  flow  was  observed  to  be 
in  the  slug  regime  [24].     It  is  assumed,  that  in  the  quality  range  0.0  to  0.1 
and  0.9  to  1.0,  the  heat  transfer  coefficient  changes  linearly  from  a  two-phase 
flow  value  to  a  single-phase  flow  value  and  is  calculated  using  linear 
interpolation  between  values  obtained  from  equations  (60)  and  (61). 


Two -Phase  Forced  Convection  with  Evaporation 


Refrigerant  enters  an  evaporator  from  an  expansion  device  at  quality  about  20 
percent  and  forms  an  annular  flow  instantly.     The  quality  increases  with  the 
proceeding  flow  and  the  annular  flow  pattern  is  maintained  until  the  quality 
reaches  about  0.90,  at  which  refrigerant  vapor  has  enough  kinetic  energy  to 
gradually  destory  the  liquid  layer  and  patches  of  dry  wall  appear. 

There  are  many  correlations  for  calculating  the  forced  convection  evaporation 
heat  transfer  coefficient  (nucleation  effect  for  mass  fluxes  taking  place  in 
heat  pump  evaporators  at  qualities  greater  than  20  percent  is  insignificant, 
even  if  it  happens).    Anderson,  Rich,  and  Geary  [26]  performed  laboratory  mea- 
surements to  obtain  the  forced  convection  evaporation  heat  transfer  coefficient 
under  conditions  closely  approximating  those  existing  in  refrigerant  evaporators. 
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They  also  examined  six  generally  known  correlations  checking  heat  transfer 
coefficient  predictions  with  their  laboratory  results.     They  found  best  agree- 
ment with  Pierre's  [27]  correlation,  which  correlate  reasonably  well  with 
experimental  data  at  high  mass  velocities  (above  10^  lb/(h  •  ft^))  and  provided 
an  excellent  agremeent  at  lower  mass  velocities. 

Pierre's  correlation,  based  on  experimental  work  with  refrigerants  R-12  and 
R-22,  has  the  following  form: 

Nu  =V-^  =  0-0009  .      G  '  D      .    r  '        '  r  ^  (66) 


where  D  =  inner  diameter  of  the  tube 

G  =  refrigerant  mass  flux 
h  =  heat  transfer  coefficient 
ifg  =  latent  heat  of  evaporation 

J  =  the  mechanical  equivalent  of  heat 
kL  =  thermal  conductivity  of  liquid  refrigerant 

L  =  length  of  the  tube 
Ax  =  refrigerant  vapor  quality  change  in  the  tube 
Ml  ~  absolute  viscosity  of  liquid  refrigerant 

Pierre's  correlation  is  applicable  for  a  Nusselt  number  value  within  the 
limit: 

28.5  <  Nu  <  753  (67) 

which  falls  in  the  range  of  heat  pump  evaporator  operation,  and  for  a  quality 
range  from  0.15  to  0.90.     For  refrigerant  quality  greater  than  0.09,  similar  to 
the  condensation  case,  the  evaporative  heat  transfer  coefficient  may  be  obtained 
by  averaging  between  the  values  yielded  by  equation  (66)  and  the  single-phase 
forced  convection  correlation,  equation  (60). 

5.4.6     Forced  Convection  Heat  Transfer  at  the  Air-Side  of  a  Flat-Finned  Tube 

In  order  to  evaluate  the  forced  convection  heat  transfer  outside  a  flat-finned 
tube  (term  4  of  equation  (59)),  the  total  exterior  surface  area,  Aq ,  the  fin 
area,  Af,  the  air-side  heat  transfer  coefficient,  h(,^Q,  and  fin  efficiency,  ((>, 
have  to  be  known. 

From  the  number  of  air-side  heat  transfer  correlation  available  in  the  literature, 
the  one  proposed  by  Briggs  and  Young  [28]   is  most  applicable  here.     This  corre- 
lation was  developed  after  extensive  tests  on  18  tube  banks  of  different  fin 
geometry.     A  regression  analysis  of  the  test  data  yielded  the  following  equation: 


h,  ^  •  0.681  0.333  .     ,  _  .  0.2        /„  ^  0.1134 

Nu=-^^^^-  ^^o.l34.Rea  -Pra  (f)        *    (f  '  ^^^^ 

where        Dq  =  outside  tube  diameter 


^max  ~         mass  flux  at  minimum  cross  section 

hj,  Q  =  air-side  mean  convective  heat  transfer  coefficient 
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=  air  thermal  conductivity 


y     •  Cp 

Pro  =  —  —  ,  Prandtl  number 

'^a 

G        •  D 

Rea  =   —  ,  Reynolds  number 

^a 

t  =  fin  thickness 
y  =  fin  height 

z  =  distance  between  adjacent  fins 


The  geometry  parameters  affecting  the  heat  tranfer  are  illustrated  in  figure 
22.  The  adequacy  of  equation  (68)  was  further  verified  by  Jones  and  Russell 
[29]. 

The  addition  of  fins  to  the  tubes  greatly  increases  the  outer  heat  transfer 
area  but  at  the  expense  of  decreasing  the  mean  temperature  difference  between 
the  surface  and  the  air  stream.     The  parameter,  fin  efficiency,  (() ,  is  used  to 
rate  the  thermal  effectiveness  of  a  fin.     As  mentioned  in  section  5.4.1,  it  is 
assumed  in  this  study  for  heat  transfer  analysis,  that  each  tube  is  served  by  a 
circular-plate  fin  of  equivalent  surface  area,  as  per  figure  21.     Gardner  [30] 
solved  the  diferential  equation  for  describing  the  temperature  distribution  in 
a  circular  fin  and  presented  fin  efficiency  curves  in  terms  of  two  parameters. 

/2  •  h^  ^  \  0.5 

Exact  theory  results  are  correlated  well  by  the  following  equation  proposed  in 
[3]  (see  also  figure  23): 


8  D  j5      2        2  .  h 

4.=    Z     (Ai,i  +A2,i  -2-   +  A3,i(-^)    )(y(-__^)        )  (69) 
i=l  ^t  ^t  ^f  ^ 

where  q  =  air-side  convective  heat  transfer  coefficient 

kf  =  fin  material  thermal  conductivity 

The  geometry  parameters  are  indicated  in  figure  22.     The  coefficients,  A^^i, 
are  given  in  table  1. 

5.4.7    Overall  Heat  Transfer  Coefficient  for  a  Wet  Finned  Tube 


Wet  finned  tube  analysis  is  applicable  to  an  evaporator  when  its  temperature  is 
below  the  dew  point  temperature  of  ambient  air.     In  such  a  case,  moisture  is 
being  removed  from  an  air  stream  and  is  transferred  to  an  evaporator  external 
surface.     If  the  evaporator  temperature  is  above  32°F,  a  water  film  flows  down 
the  fin  under  force  of  gravity.     If  the  exterior  evaporator  temperature  is 
below  32°F,  frost  is  accumulated. 
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Do  =  TUBE  OUTSIDE  DIAMETER 
Dt  =  FIN  TIP  DIAMETER 
y   =  FIN  MIGHT 
t    -  FIN  THICKNESS 

z    =  DISTANCE  BETWEEN  ADJACENT  FINS 


Figure  22.     Cross  section  of  a  flat-finned  tube  indicating  parameters 
affect  the  air-side  heat  transfer  coefficient 
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Figure  23.    Efficiency  for  a  circular-plate  fin  of  uniform  thickness. 

Comparison  of  exact  theory  results  with  those  obtained 
by  equation  (69) 
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Table  1.    Coefficients  to  be  Used  in  Correlation  for  Fin  Efficiency 
(Equation  (69)) 


i 

^l,i 

A,2i 

^3,1 

1 

1.0 

0.0 

0.0 

2 

-0.22920E-01 

-0.13755E+00 

0.20130E- 

■01 

3 

0.16106E+00 

0  81890E- 

•01 

-0.11440E- 

■01 

4 

-0.64975E-I-O0 

-0.55500E- 

■01 

-0.28753E- 

•01 

5 

0.53491E+00 

0.18040E- 

■01 

0.42477E- 

•01 

6 

-0.1928 6 E+OO 

0.36494E- 

•03 

-0.20335E- 

•01 

7 

0.32564E-01 

-0.10660E- 

■02 

0.40947E- 

■02 

8 

-0.20972E-02 

0.12410E- 

-03 

-0.29673E- 

■03 

65 


Looking  at  the  heat  transfer  phenomena  from  refrigerant  to  air,  one  can  realize 
that  the  dehumldif ication  process  will  alter  the  heat  transfer  situation  on  the 
external  surface  of  the  finned  tube,  while  other  processes  in  the  tube  and 
refrigerant  stay  unaffected  and  are  governed  by  relations  already  proposed  in 
previous  sections.     Thus  only  processes  connected  with  dehumidif ication  need  to 
be  discussed. 

The  heat  transfer  rate  between  the  air  stream  and  the  water  surface  is  described 
by  the  following  equation: 

dQ  =  (hc,o(Ta  -  Tw)  +  tiD.oC^a  "  Ww)ifg,w)dAo  (70) 

The  first  term  in  equation  (70)  accounts  or  sensible  heat  and  the  second  term 
accounts  for  latent  heat  transfer.     Noting  that  for  air  at  atmospheric  pressure 
the  Lewis  number, 

K  « 

Le  =  Sut2-  ,  (71) 

^D,oCpa 

is  close  to  1,  equation  (70)  assumes  the  following  form  for  a  flat  finned  tube: 

dQ  =  he,o(l  +  ^pg'^r^  "V  )(1  (1  -  <f))(Ta  "  Tw)dAo  (72) 

Cp^(T^  -  T^) 

Symbols  used  in  equations  (70),  (71),  and  (72)  denote: 

Af  =  fin  area 

Aq  =  total  external  area 

CPa  ~  specific  heat  of  air 

h(,  Q  =  air-side  forced  convection  heat  transfer  coefficient 

^D,o  ~  air-side  mass  transfer  coefficient 

^fg,w  ~  latent  heat  of  condensation  for  water 

T^  =  temperature  of  air 

=  temperature  of  liquid  water  (frost) 

Q  =  heat  transfer  rate 

W3  =  humidity  ratio  of  air 

w^,  =  humidity  ratio  of  saturated  air  at  T^  temperature 

T.       -  T 
^  =    w,ni  ^ — a  ^  efficiency 

^w,b  ^a 

^w,m  ~  insai^  temperature  of  water  film  (frost) 

^w,b  ~  temperture  of  water  film  (frost)  at  fin  base 

The  one-dimensional  heat  conduction  across  the  condensate  (frost)  film  can  be 
expressed  by  equation: 

dQ  =  hL  .  ATl  •  dAo  (73) 
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where 


_  K 

  ,  heat  transfer  coefficient  for  condensate  (frost)  film 

6 

=  thermal  conductivity  of  water  (frost) 
6  =  thickness  of  condensate  (frost)  film 
ATl  =  temperature  difference  across  condensate  (frost)  film 


Using  equations  (72)  and  (73)  and  referring  to  equation  (59),  the  following 
relation  for  overall  heat  transfer  coefficient  for  a  wet  finned  tube  can  be 
derived: 

-1 


U  = 


hiAp,i 


^d,i^p,i 


A  k 
"^p  ,rorp 


(,^iffi,w("a-"»))(,.Af(i-  *)) 


CO 


Cp  (T  -  T  ) 
^a^  a 
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where  symbols  used  are  defined  as  in  equations  (59)  and  (72). 

In  the  above  formulation  of  the  overall  heat  transfer  coefficient,  it  is  assumed 
that  the  temperature  difference  across  the  liquid  film  (frost)  is  uniform.  The 
calculated  value  for  wet  fin  efficiency  is  affected  by  the  change  of  the  air- 
side  heat  transfer  coefficient  caused  by  air  passage  geometry  alternation  by 
liquid  (frost)  accumulation,  and  by  released  latent  heat  of  condensation  (sub- 
limation).    The  effect  of  water  (frost)  conductance  on  fin  efficiency  is 
neglected.     In  summary  the  heat  transfer  phenomena  that  occurs  during 
dehumidif ication  on  the  air-side  may  be  itemized  as  follows: 

(1)  layer  of  wet  (frost)  offers  additional  heat  flow  resistance 

(2)  air-side  heat  transfer  resistance  is  decreased  due  to  effect  of 

condensation 

(3)  air-side  heat  transfer  coefficient,  h^,  q,  since  it  is  sensitive 

to  external  surface  geometry,  has  an  increased  value  (see 
equation  (68)) 

(4)  fin  efficiency  has  decreased  value  as  h^,  q  is  increased  (see 

figure  23) 

(5)  the  cross  sectional  area  of  the  air  flow  passage  between  the  fins 

has  decreased,  decreasing  the  flow  rate  and  the  heat  exchange. 

In  order  to  evaluate  water  (frost)  layer  thickness,  consider  the  mass  transfer 
equation: 

^a,d  •  dwa  =  -  hD^o(Wa  -  w^)dAo  (75) 

For  the  Lewis  number  equal  to  1  and  for  a  finned  tube  differential  surface  area, 
equation  (75)  assumes  the  following  form: 

^Pa  ^o 
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Variation  of  air  humidity  ratio  can  be  calculated  by  integrating  equation  (76), 
which  yield&t 


-  ^c,o  '  -^(1  - 

Wa,e  =  Wa^i  -  (w^^i  -  wj(l  -  exp   Cv    '  m   ^  ^^^^ 

"a        a  ,d 

The  rate  of  moisture  removal  per  unit  area,  R,  can  be  now  calculated: 

R  =  ma^d(wa,i  "  ^a,e^/^  (78) 

where         d  ~  '^^ss  flow  rate  of  dry  air 

Wg  e  =  humidity  ratio  of  air  at  tube  row  exit 
Wg^i  =  humidity  ratio  of  air  at  tube  row  inlet 

If  the  evaporator  temperature  is  below  the  freezing  point,  moisture  removed 
from  the  air  stream  accumulates  on  the  evaporator  external  surface  in  the  form 
of  frost.     Its  thickness,  6f ,  can  be  evaluated  by  integrating  the  rate  of 
moisture  removal  with  respect  to  time,  i.e.: 


t 

6f  =  /      1    dt  (79) 
0  Pf 

where         t  =  time 

R  =  rate  of  moisture  removal  per  unit  area 
6f  =  frost  layer  thickness 
=  frost  density 

In  case  of  evaporator  temperature  above  32° F,  condensate  flows  down  on  the  fin. 
Assuming  no  drag  on  the  liquid  layer,  the  condensate  thickness,  6^,  can  be 
found  by  solving  the  Navier-Stokes  equation.     The  resulting  expression  for  a 
finned  tube  situation  is: 


y    •  R' 
6f  =  1.083  (-^  -) 

g  •  P 


(80) 


where         g  =  gravitational  acceleration 

R'  =  condensation  rate  per  unit  width  of  a  fin 
yjjj  =  water  dynamic  viscosity 
=  water  density 

5.4.8    Pressure  Drop  in  a  Tube 

As  expressed  by  equation  (25) j  the  total  pressure  drop  experienced  by  a  flowing 
fluid  results  from  pressure  drops  due  to  friction,  momentum  change,  and  gravity. 
In  an  actual  heat  pump  heat  exchanger,  pressure  drop  due  to  gravity  effect  is 
very  small  and  may  be  neglected.     Only  pressure  drop  due  to  friction  and  due  to 
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momentum  change  will  be  considered  with  distinction  between  different  flow 
patterns  in  a  tube. 


Single-Phase  Flow 

Frictional  pressure  drop  for  a  single-phase  in  a  tube  can  be  calculated  by  the 
Fanning  equation  with  the  Fanning  friction  factor  as  per  equations  (43)  and  (44): 

AP  "        •        •  ^  (81) 
D  •  p 

f  =  0.046  Re"0-2  (82) 

Pressure  drop  due  to  momentum  change  can  be  calculated  by  the  following  equation; 

dP  =  _  g2  dv  (83) 
dL  dL 

where  G  =  refrigerant  mass  flux 

L  =  tube  length 

V  =  refrigerant  specific  volume 
Two-Phase  Flow  with  Condensation 

Frictional  pressure  drop  for  two-phase  flow  with  condensation  can  be  calculated 
by  the  method  proposed  by  Lockhard  and  Martinelli  [25].     They  performed  a  semi- 
empirical  study  of  adiabatic  two-phase  flow  with  air  and  different  liquids 
including  benzene,  kerosene,  water,  and  various  oils  in  pipes  varying  in  diam- 
eter form  0.586  to  1.017  inch.     They  related  pressure  drop  of  two-phase  flow 
to  the  pressure  drop  of  the  liquid  portion  of  the  flow  flowing  alone  in  the 
pipe,  by  the  dimensionless  parameter  X^-|-,  i.e.: 

=  f(Xtt)  =  *  APtp  =  APl  •  *  (8^) 


where      APl  =  frictional  pressure  drop  of  the  liquid  portion  of  two-phase  flow 
flowing  alone  in  the  tube 
AP^p  =  frictional  pressure  drop  of  two-phase  flow 
9  =  correction  factor  for  two-phase  pressure  drop 
Xj.^.  =  as  given  by  equation  (64) 

The  pressure  drop,  APl,  can  be  calculated  by  the  single-phase  pressure  drop 
relation  with  calculated  liquid  Reynolds  number  and  friction  factor  as  follows: 

Re^„  T  =  <^  -  x)G  »  D  (85) 

f,    ,  =  0.046  •  Re,      "0-2  (86) 
tp,L  tp,L 

where  ftp,L  =  friction  factor  for  the  liquid  portion  of  two-phase  flow  flowing 
alone  in  the  pipe 
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^^tp,L  ~  Reynolds  number  for  the  two-phase  liquid  portion  flowing  alone 
in  the  pipe 
X  =  quality 
PL  ~  liquid  dynamic  viscosity 

A  correction  factor  for  two  phase  pressure  drop,  0,  can  be  correlated  by  the 
following  equation: 

Z    A^Clog  X^^)^"^ 

$  =  10  i=l  (87) 

where       ^1  ~ 

A2  =  -0.87917 
A3  =  0.14062 
A4  =  0.0010417 
A5  =  -0.00078125 


Combining  equations  (84),   (85),  (86),  and  (87),  the  two-phase  pressure  drop 
equation  assumes  the  form: 

AP^p  =  2  f^p^L  •  -  x^)L  .  $/(D  .  pl)  (88) 

The  pressure  drop  due  to  momentum  change  for  separated  two-phase  flow  can  be 
estimated  by  the  following  equation: 

2  2 

jr)           OA      VtT  •  X  Vt  (1  -  x) 

dP=-G2d     (-V   +-Jf-  _)  (89) 

dL              dL         a  Ql  -  a; 

where         x  =  quality 

vl  =  specific  volume  of  liquid 
vy  =  specific  volume  of  vapor 
a  =  void  fraction 

Void  fraction,  a,  percent  of  tube  filled  with  vapor,  was  shown  by  Lockhard  and 
Martinelli  to  be  a  function  of  X^^-  under  any  flow  conditions  for  separated  flow 
with  both  phases  turbulent.     Wallis  [31]  correlated  their  results  in  the 
following  form: 

n  «  -0,378 

a  =  (1  +  X^.^0-»)  (90) 

This  expression  was  found  to  correlate  well  with  data  presented  in  [25]  for 
values  of         _<  10.     For  X^.,.  greater  than  10,  another  curve  fitted  formula  is 
used: 

a  =  0.823  -  0.157  .  In  X^t  (91) 
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Two-Phase  Flow  with  Evaporation 


The  pressure  drop  for  two-phase  flow  with  evaporation  can  be  calculated  by  the 
correlation  proposed  by  Pierre  [32],  based  on  his  experiments  with  refrigerants 
R-12  and  R-22.     Pierre's  correlation  has  the  following  form: 

AP  =  (f  L  +  A2C)g2  .  Vjj^  (92) 

where         f  =  friction  factor  (calculated  by  equation  (93)) 
=  mean  quality 
Ax  =  quality  change 

Vpj  =  vl  +  Xnj(vY  -  vl)  ,  mean  specific  volume 

The  friction  factor  to  be  used  in  equation  (92)  was  correlated  by  Pierre  from 
his  experimental  data  by  the  following  empirical  equation: 

f  =  0.0185  (Re   h  )0-25  (93) 

J  •  ifg  •  Ax 

where        Re  =  ^  *  P  ,  Reynolds  number 

The  correlation  proposed  by  Pierre  is  in  the  conventional  format  for  the  single 
pressure  drop  formula.    The  first  term  of  equation  (92)  is  for  frictional  pres- 
sure drop  while  the  second  is  for  pressure  drop  due  to  change  of  momentum.  The 
formula  for  the  friction  factor  contains  conventionally  the  Reynolds  number 
divided  by  the  term  J  •  ifg  •  Ax/L,  referrred  to  as  a  boiling  number,  making 
friction  factor  sensitive  to  vapor  generation  rate  at  the  vapor-liquid  inter- 
face.    Pierre's  correlation  was  verified  in  [26]  providing  overall  better 
agreement  with  experimental  data  for  pressure  drop  of  refrigerants  R-12  and 
R-22  than  the  other  popular  correlation  of  Martinelli  and  Nelson  [33] . 

5.5    MODELING  OF  ADDITIONAL  HEAT  PUMP  COMPONENTS 

In  the  previous  sections,  modeling  of  the  main  heat  pump  components  has  been 
discussed.     The  analysis  included  a  compressor,  a  constant  flow  area  expansion 
device,  a  condenser,  and  an  evaporator.     Additional  components  that  have  to  be 
considered  for  a  more  accurate  heat  pump  performance  prediction  by  the  model 
are:     vapor  suction  and  discharge  lines,  liquid  line,  and  a  four-way  valve. 
Their  modeling  is  briefly  explained  below. 

Vapor  Lines 

A  suction  line  connects  an  evaporator  with  a  compressor.     A  compressor  and  a 
condenser  are  connected  by  a  discharge  line.     In  both  lines,  single-phase  flow 
exists  in  the  form  of  either  saturated  or  superheated  vapor.     In  heat  pump 
systems  vapor  lines  are  usually  insulated.     Heat  transfer  rates  between  refri- 
gerant vapor  flowing  in  these  pipes  and  ambient  air  can  be  calculated  by  a 
general  heat  transfer  equation  for  a  circular  duct  with  insulation  [10]. 
Single-phase  forced  convection  is  assumed  inside  the  tube  (equation  (60))  and 


71 


free  convection  is  assumed  outside  the  tube.     The  following  equation  is  used 
for  calculation  of  the  free  convection  heat  transfer  coefficient  for  a 
horizontal  pipe  [38] : 

A  rri     0  •  2  5 

h  =  0.27  (^)  (94) 

for  Grashof  numbers  from  10^  to  109 

where       AT  =  temperature  between  tube  wall  and  air 
Dq  =  tube  outside  diameter 

Pressure  drop  in  vapor  lines  can  be  calculated  by  the  single-phase  pressure 
drop  equation  (43). 

Liquid  Line 

A  liquid  line  connects  an  evaporator  and  a  condenser.     This  line  is  filled  with 
subcooled  liquid  or  low  quality  two-phase  flow.     The  pressure  drop  can  be  cal- 
culated by  equations  (43)  or  (88),  depending  on  the  flow  pattern.     Heat  loss 
from  the  liquid  line  to  the  ambient  air  is  neglected. 

Four-Way  Valve 

The  main  function  of  a  four-way  valve  is  to  direct  refrigerant  flow  from  the 
compressor  to  the  indoor  or  outdoor  coil  depending  on  the  mode  of  operation 
(heating  or  cooling).     The  side  effects  of  flow  through  a  four-way  valve  are 
changes  in  the  refrigerant  thermodynamic  state  due  to  pressure  drop  and  heat 
exchange.     Assuming  an  adiabatic  exterior  wall,  all  heat  lost  by  the  discharge 
refrigerant  is  gained  by  the  suction  refrigerant.     Heat  transfer  rate  and  refri- 
gerant pressure  drop  in  the  valve  can  be  evaluated  by  formulas  similar  to  equa- 
tions (20)  and  (30),  respectively.     Heat  transfer  and  pressure  drop  parameters 
in  these  formulas  have  to  be  found  empirically. 

5.6     REFRIGERANT  MASS  INVENTORY  IN  A  HEAT  PUMP 

The  mass,  M,  of  a  substance  occupying  a  known  volume,  V,  may  be  determined  by: 

M  =  /  p  •  d¥  (95) 
¥ 

where  p  =  local  density 

In  reality,  the  mass  of  refrigerant  in  the  system  can  be  fund  by  estimating 
masses  of  the  refrigerant  in  each  system  component  and  adding  them  up.  For 
this  purpose,  equation  (95)  can  be  written  in  the  form: 

M  =  E  Mi  (96) 

Mi  =  Vi  .  p^^i 
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where        Mj^  =  mass  of  refrigerant  in  particular  component  1 
Vj^    =  internal  volume  of  component  i 
Pm,i  ~  mean  fluid  density  in  component  i 

The  heat  pump  components  taken  into  account  in  the  mass  inventory  calculations 
are  shovm  in  figure  24.     The  refrigerant  phase  in  each  of  the  components  as 
indicated  in  the  figure  is  based  on  the  following  consideration: 

Discharge  line  -  receives  and  is  filled  with  superheated  vapor  from  the 
compressor . 

Condenser  -  receives  superheated  vapor  from  a  discharge  line.     In  the 
course  of  passage  through  the  condenser  tubes,  vapor  temperature  is  brought 
to  saturation  temperature.     Starting  at  this  point,  a  thin  condensed  liquid 
layer  forms  on  the  tube  walls.     Depending  on  the  mass  flux,  this  liquid 
film  may  be  swept  and  entrained  within  the  vapor  as  a  mist  forming  a  dis- 
persed flow.    With  more  condensed  vapor,  the  velocity  decreases,  a  perman- 
ent liquid  film  is  formed  on  the  wall,  and  flow  proceeds  further  as  annular 
flow  or  semiannular  flow.     The  quality  decreases  and  flow  slows  down  in 
the  direction  of  flow.     Gravity  forces  result  in  a  stratifying  effect  and 
liquid  flowing  in  the  bottom  of  the  tube  periodically  wets  the  upper  tube 
wall  (slug  flow).     If  full  condensation  is  reached,  there  will  be  subcooled 
liquid  at  the  condenser  exit. 

Thus  there  are  four  two-phase  flow  patterns  in  the  refrigerant  condensers, 
i.e.: 

dispersed  flow,  annular  flow,  semi-annular  flow,  slug  flow 

as  identified  by  Traviss  and  Rohsenow  [34].     Soliman  and  Azer  [35] 
identified  nine  flow  regimes;  however,  five  of  them  may  be  considered  as 
transition  regimes  and  can  be  put  into  the  above  basic  four  categories. 
The  effect  of  return  bends  connecting  condenser  tubes  was  investigated 
in  [36].     It  was  found  that  this  effect  is  insignificant  and  that  the 
annular  flow  pattern  "recovers"  almost  immediately  in  a  tube  after  passing 
a  return  bend. 

Liquid  line  -  receives  refrigerant  from  a  condenser  and  delivers  it  to  an 
expansion  device.     Entering  the  small  bore  restrictor  at  the  end  of  the 
liquid  line,  the  flow  experiences  a  sudden  contraction.     The  entrance  pre- 
sure  drop  depends  upon  the  mass  flow  rate  and  fluid  density.     This  along 
with  capillary  sensitivity  to  the  amount  of  subcooling,  provides  a  flow 
controlling  mechanism  for  the  system.     It  is  common  practice  to  design  a 
capillary  tube  so  the  refrigerant  mass  flow  rate  through  it  would  balance 
with  the  mass  flow  rate  through  the  compressor  when  a  liquid  seal  exists 
before  the  capillary  tube.     Temporary  existence  of  vapor  before  the  capil- 
lary tube  would  decrease  the  mass  flow  rate  and  build  up  the  pressure. 
This  would  provide  better  conditions  for  condensation  in  the  condenser  and 
thus  cause  a  return  to  a  subcooled  liquid  or  very  low  quality  two-phase 
flow  in  the  liquid  line. 
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LIQUID  LINE 


LIQUID  PHASE  OR  LOW 
QUALITY  TWO-PHASE 


Figure  24.    Refrigerant  phase  in  heat  pump  components 
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Capillary  tube  -  due  to  small  internal  volume  may  be  disregarded  for  mass 
inventory  purposes. 

Evaporator  -  receives  refrigerant  from  a  capillary  tube  in  the  form  of 
homogeneous  flow  of  about  20  percent  quality,  at  sonic  or  close  to  sonic 
velocity.     This  flow  experiences  a  sudden  deceleration  and  is  assumed  to 
form  an  annular  or  semi-annular  flow  pattern  very  quickly  having  a  void 
fraction  about  0.8  [37] .     The  liquid  phase  flows  along  the  tube  in  the 
form  of  an  annulus  and  the  vapor  phase  flows  in  the  core.     This  type  of 
flow  prevails  in  most  of  the  evaporator.     With  increasing  quality,  the 
velocity  of  the  vapor  increases  and  becomes  eventually  high  enough  to 
entrain  small  liquid  drops  and  tear  the  liquid  film  apart.     This  results 
in  a  dry  mist  flow  at  qualities  close  to  1.     Finally,  the  flow  may  reach 
the  evaporator  exit  in  a  superheated  vapor  state. 

Suction  line  -  receives  refrigerant  from  the  evaporator  as  a  mist  flow  or 
superheated  vapor  flow  and  the  same  flow  pattern  exists  at  the  exit. 

Accumulator  -  receives  the  refrigerant  from  a  suction  line  after  some 
heat  was  added  in  the  four-way  valve.  Saturated  vapor  or  superheated 
vapor  may  be  expected.  (This  does  not  account  for  refrigerant  liquid 
that  may  be  accumulated  on  the  accumulator  bottom.) 

Compressor  can  -  receives  refrigerant  from  the  accumulator  as  a  saturated 
or  superheated  vapor  which  is  further  superheated  due  to  heat  transfer 
from  an  electric  motor,  discharge  line,  and  compressor  cylinder  body. 

According  to  the  above  analysis,  refrigerant  in  the  accumulator,  compressor 
can,  suction  line,  and  discharge  line  is  in  (or  close  to)  a  homogeneous  vapor 
state.     In  the  liquid  line,  homogeneous  liquid  flow  or  two-phase  flow  exists 
but  in  a  state  so  close  to  saturation  that  both  liquid  and  vapor  would  travel 
at  about  the  same  velocity.     For  mass  inventory  purposes,  the  flow  may  be 
assumed  to  be  homogenous.     The  mass  of  refrigerant  in  these  components  can  be 
then  calculated  in  a  similar  manner  since  refrigerant  densities  in  the  compo- 
nents described  above  will  be  known  as  a  result  of  the  simulation  program 
iteration  process. 

In  the  case  of  both  coils,  the  refrigerant  flow  in  part  is  homogeneous;  however, 
in  most  of  the  coil  some  type  of  annular  flow  prevails.     For  a  separated  flow 
regime  like  an  annular  flow,  density  of  the  flowing  fluid  results  from  vapor 
and  liquid  densities  and  the  fractions  of  tube  volume  occupied  by  the  liquid 
and  vapor  phase: 

Pm  =  „^  =  —  ^^r^   =  Pl(1  -  a)  +  Pv  •  «  (97) 

t  ^t 

where         M  =  mass 

V  =  occupied  volume 
p  =  density 

a  =  Vv/Vt,  void  fraction 
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subscripts  refer  to: 


L  =  liquid  phase 
m  =  mean  value 
t  =  total  value 
V  =  vapor  phase 


It  should  be  pointed  out  that  for  a  two-phase  flow  of  given  quality,  mean 
density  depends  not  only  on  the  thermodynamic  parameters  affecting  densities  of 
the  liquid  and  vapor,  but  also  on  the  ratio  of  mean  velocities  of  the  vapor 
and  liquid  referred  often  as  the  slip  ratio.     This  can  easily  be  noticed  by 
looking  at  the  equation  for  void  fraction,  which  can  be  derived  considering 
mass  flow  rates  of  each  phase: 


a  =   i   (98) 

1  -  X  Vy  py 
^       Vl  PL 

where         V  =  velocity 


Calculation  of  the  void  fraction  has  received  much  less  attention  than  the 
calculation  of  heat  transfer  or  pressure  drop,  although  all  three  of  these 
phenomena  are  undoubtedly  inter-related.     The  most  often  referred  method  for 
void  fraction  calculation  is  that  of  Lockhard  and  Martinelli  [25]  based  on 
their  experimental  void  fraction  data  reported  along  with  pressure  drop  results. 
Lockhard  and  Martinelli  correlated  void  fraction  with  the  dimensionless  param- 
eter, Xj-t  •     Their  experiment  dealt  with  adiabatic  flow;  however,  in  the  conclu- 
sions of  their  paper  they  suggested  that  their  pressure  drop  correlation  could 
be  used  for  prediction  of  pressure  drop  during  evaporation  and  condensation  as 
well.     On  similar  grounds,  their  void  fraction  data  should  be  applicable  beyond 
the  adiabatic  case. 


There  is  another,  earlier,  method  available  for  void  fraction  calculation  by 
Martinelli  and  Nelson  [33],  derived  with  the  same  assumptions  for  water/steam 
evaporating  flow.     The  method  is  similar  to  the  Lockhard-Martinelli  method  so 
the  Lockhard  and  Martinelli  method  was  used  in  this  study  for  calculation  of 
the  void  fraction  in  both  evaporator  and  condenser  coils. 

Correlations  for  void  fraction  based  on  [25]  have  already  been  given  in  section 
5.4.8  (equations  (90)  and  (91)). 

a  =  (1  +  Xj.,.^*^)  "  ^-^^^         for  Xj.^.  <  10  (99) 

a  =  0.823  -  0.157  In  Xtt  for  Xtt  >  10  (100) 

Parameter  Xj-t  >  as  defined  by  equation  (64),  is  a  function  of  flow  quality, 
specific  volume,  and  viscosities  of  the  liquid  and  vapor,  and  is  not  sensitive 
to  the  mass  flux  which  in  turn  affects  the  slip  ratio. 
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Experimental  data  of  Staub  and  Zuber  [37]  for  evaporating  refrigerant  R-22 
indicate  that  void  fraction  increases  with  increased  mass  flow  rate.  Compari~ 
son  of  their  data  with  void  fraction  predicted  by  the  Martinelli-Nelson  method 
showed  a  discrepancy  which  may  lead  to  underestimation  of  the  two-phase  flow 
mean  density  by  as  much  as  300  percent.     (Note  that  a  small  difference  in  void 
fraction  results  in  a  large  error  in  mean  density  prediction  as  density  of  the 
vapor  and  liquid  are  of  different  orders  of  magnitude.)     The  results  of  Staub 
and  Zuber  also  indicated  that  approximate  agreement  with  the  Martinelli-Nelson 
method  could  be  obtained  at  high  mass  fluxes,  more  than  five  times  higher  than 
those  observed  in  heat  pump  evaporators. 

In  discussing  the  accuracy  of  a  refrigerant  mass  prediction  method,  it  should 
be  noted  that  it  also  depends  upon  an  accurate  measurement  of  the  internal 
volume  of  the  system.     Internal  volume  can  be  easy  calculated  for  straight 
pipes.     For  valves,  bends,  etc.,  it  can  be  only  approximated. 

Mass  inventory  of  a  heat  pump  system,  as  explained  in  section  5.1.2,  is  intended 
to  be  used  for  iteration  of  refrigerant  superheat  at  the  compressor  can  inlet. 
The  discussions  above  indicated  that  mass  inventory  calculations  may  have  errors 
from  two  different  sources.     These  sources  are:     inaccurate  density  prediction 
in  the  two-phase  region,  and  inaccuracte  internal  volume  measurement.  However, 
these  inaccuracies  are  not  of  great  importance  and  can  be  ignored,  since  the 
program  does  not  really  need  prediction  of  the  absolute  value  of  refrigerant 
mass  in  the  system,  but  rather  requires  sensitivity  in  relative  mass  predictions 
in  a  heat  pump  with  change  of  operating  conditions .     This  requirement  should  be 
satisfied  by  the  Lockhard-Martinelli  method  and  very  precise  knowledge  of  the 
internal  volume  should  not  be  that  important.     Since  refrigerant  vapor  superheat 
at  the  compressor  can  inlet  is  known  for  a  heat  pump  operating  in  the  cooling 
mode  at  design  conditions  (outdoor  temperature  90°F,  indoor  dry  bulb  temperature 
SO'F,  and  wet  bulb  67°F)  and  does  not  have  to  be  iterated  upon,  refrigerant 
mass  in  the  system  can  be  calculated  at  this  operating  condition.     This  calculated 
mass  of  refrigerant  may  then  be  used  as  input  data  for  refrigerant  vapor  super- 
heat iteration  at  other  outdoor  or  indoor  air  conditions  at  which  heat  pump 
simulation  results  are  required. 
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6.     MODEL  VERIFICATION 


A  verification  of  the  computer  model  has  been  done  by  comparison  of  computer 
predictions  with  performance  data  obtained  from  laboratory  tests  conducted  in 
NBS  environmental  chambers .     Four  heat  pumps  were  available  for  computer  model 
evaluation.     They  are  denoted  here-ln  as  systems  1,  2,  3,  and  4.     System  1, 
system  2,  and  system  3  consisted  of  the  same  compressor,  outdoor  coil,  con- 
necting tubing,  and  expansion  device.     The  only  variable  element  was  an  indoor 
section  (containing  an  Indoor  fan  and  coil).     In  the  case  of  system  2,  the 
Indoor  unit  was  an  original  match  for  the  rest  of  the  components.     For  system 
1  and  for  system  3,  the  Indoor  sections  were  undersized  and  oversized,  respec- 
tively.    System  4  was  a  matched  indoor/outdoor  system  of  another  manufacturer. 
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6.1     LABORATORY  TESTS 


An  identical  series  of  four  tests  was  conducted  in  the  cooling  mode  on  each  of 
the  four  systems.  The  series  consisted  of  a  test  at  an  outdoor  temperature  of 
105°F  and  tests  A,  B,  and  C  of  the  NBS  recommended  test  procedure  [39] .  Ther- 
mal conditions  during  the  above  tests  covered  practically  the  whole  range  of 
the  cooling  mode  operating  conditions  for  a  heat  pump  working  in  the  field. 
Additionally,  heat  pump  system  2  was  tested  in  the  heating  mode  at  temperature 
conditions  recommended  by  [40] . 

Cooling  and  heating  tests  results  are  presented  in  table  2  and  table  3, 
respectively.     Symbols  used  in  these  and  following  tables  in  this  section  are 
consistent  with  symbols  used  in  the  main  program  MAIN  (refer  to  appendix  I)  and 
have  the  following  meaning: 


COP      -  coefficient  of  performance 

ELUSE  -  energy  consumption  rate 

P3        -  refrigerant  pressure  at  compressor  inlet 

P8        -  refrigerant  pressure  at  compressor  outlet 

RHOA    -  relative  humidity  of  outdoor  air 

RHRA    -  relative  humidity  of  indoor  air 

RMASS  -  refrigerant  mass  flow  rate 

TOA      -  temperature  of  outdoor  air 

TRA      -  temperature  of  indoor  air 

TSUP3  -  superheat  at  compressor  inlet 

QLOAD  -  heat  pump  capacity 


The  capacity  reported  in  table  2  is  the  one  measured  by  the  refrigerant  enthalpy 
method  [41].     In  this  method,  the  unit  capacity  is  derived  from  the  refrigerant 
enthalpy  change  and  refrigerant  volumetric  flow  rate  measurements  with  adjust- 
ment for  heat  added  to  the  indoor  section  by  the  indoor  fan.     Capacity  obtained 
by  the  refrigerant  enthalpy  method  is  directly  comparable  to  the  capacity 
obtained  from  the  vapor  compression  cycle  simulation. 

During  heating  mode  tests,  no  liquid  subcooling  existed  in  the  liquid  line  of 
the  heat  pump.     Vapor  bubbles  existing  in  the  flow  made  it  impossible  to  mea- 
sure volumetric  flow  rate  of  refrigerant.     In  this  case,  system  capacity  was 
determined  by  the  air-side  enthalpy  difference  method  [41] . 


6.2     COMPUTER  SIMULATION  RUNS 


Each  laboratory  test  performed  was  simulated  using  the  computer  program  HPSIM. 
Since  laboratory  tests  were  performed  first,  the  computer  simulations  could 
be  made  with  the  same  thermal  operating  conditions  that  actually  existed  in  the 
laboratory  during  the  tests. 

The  first  run  executed  for  each  system  was  at  design  operating  conditions  (test 
2)  with  refrigerant  superheat  at  the  compressor  inlet  equal  to  that  of  laboratory 
test.     The  refrigerant  charge  calculated  during  this  run  was  used  later  as  input 
data  for  other  runs  requiring  superheat  iteration. 
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Table  2.     Laboratory  Test  Results  of  Heat  Pumps  in  the  Cooling  Mode 


TRA 
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ELUSE 
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RMASS 

PJ 

TbUP-J 

no 

Po 

_ 

_ 

_ 

Btu/h 

Watt 

_ 

Ib/h 
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Z  J  .4 
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1 

1.4 

82.1 

79.5 

.19 

33203 

4511 

2.157 

534 

73.0 

1.0 

210.5 

2 

2.1 

104.9 

80.0 

.52 

34841 

5212 

1.959 

589 

88.9 

2.3 

292  .2 

2 

2.2 

94.9 

80.0 

.52 

37548 

4938 

2.228 

578 

85.0 

13.8 

259  .5 

2 

2.3 

82.1 

80.0 

.52 

40105 

4629 

2.538 

555 

79.3 

30.9 

220.2 

2 

2.4 

82.1 

79.8 

.20 

37156 

4583 

2.375 

550 

76.3 

9.4 

216.6 

3 

3.1 

105.2 

79.9 

.51 

37963 

5542 

2.007 

640 

95.5 

1.2 

305.2 

3 

3.2 

95.2 

79.8 

.52 

40946 

5276 

2.274 

622 

90.9 

12.8 

270.3 

3 

3.3 

81.9 

79.9 

.51 

42842 

4930 

2.546 

602 

85.1 

21.9 

230.8 

3 

3.4 

82.2 

79.9 

.21 

38551 

4867 

2.321 

591 

80.2 

1.1 

224.9 

4 

4.1 

104.9 

80.0 

.51 

31595 

4613 

2.007 

499 

89.3 

20.0 

316.7 

4 

4.2 

95.0 

80.0 

.51 

32640 

4351 

2.198 

483 

84.8 

25.0 

280.1 

4 

4.3 

82.2 

80.0 

.51 

33656 

3982 

2.476 

449 

76.6 

32.5 

236.2 

4 

4.4 

82.1 

80.0 

.18 

31437 

3971 

2.320 

445 

75.2 

25.9 

233.6 
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Table  3.    Laboratory  Test  Results  of  System  2  in  the  Heating  Mode 


SYSTEM 

TEST 

TOA 

RHOA 

TRA 

QLOAD 

ELUSE 

COP 

P3 

TSUP3 

P8 

°F 

°F 

Btu/h 

Watt 

psia 

°F 

psia 

2 

2.5 

46.6 

.73 

70.2 

42648 

4763 

2.624 

71.4 

0.0 

239.9 

2 

2.6 

35.0 

.82 

69.9 

35553 

4471 

2.330 

60.4 

0.0 

218.2 

2 

2.7 

17.1 

.70 

70.0 

24958 

3956 

1.848 

45.6 

0.0 

187.9 
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Cooling  simulation  results  are  presented  in  table  4.     Heating  simulation 
results  are  presented  in  table  5.     A  comparison  of  results  is  given  in  table  6. 
This  table  shows  that  the  maximum  discrepancy  found  between  laboratory  test 
results  and  computer  simulation  results  was: 

capacity  -  3.4  percent 

energy  consumption  rate  -  7.3  percent 

COP  -  5.5  percent 

The  largest  deviation  of  refrigerant  conditions  was: 

refrigerant  mass  flow  -  5.2  percent 
refrigerant  suction  pressure  -  3.8  psi 
refrigerant  superheat  at  suction  -  8.4°F 
refrigerant  discharge  pressure  -  23  psi 

Agreement  between  results  was  generally  very  good.     The  only  serious  difference 
was  in  the  case  of  refrigerant  discharge  pressure  where  the  discrepancy  was 
found  to  be  as  much  as  23  psi,  which  corresponds  to  about  7.5°F  of  difference 
in  refrigerant  saturation  temperature  in  the  condenser.     This  disagreement  took 
place,  however,  only  in  the  heating  mode  of  operation  when  the  indoor  coil  was 
the  condenser.     Discharge  pressure  results  for  cooling  operation  agreed  much 
better. 

In  the  course  of  the  analyses  of  this  23  psi  pressure  difference  the  following 
conclusions  were  made: 

-  there  was  a  good  agreement  between  refrigerant  suction  pressure  results 
and  capacity  results  for  cooling  mode  operation,  which  indicates  that 
the  indoor  coil  is  simulated  well  while  operating  as  an  evaporator. 

-  there  was  a  significant  difference  in  the  heat  transfer  rates  in  the 
performed  tests  in  the  indoor  and  outdoor  coil  working  as  a  condenser. 
Data  indicated  that  the  heat  flux  was  twice  as  large  in  the  indoor  coil 
than  in  the  outdoor  coil,  which  has  considerably  more  surface  area  and 
volume.     Additionally,  the  refrigerant  mass  flow  rate  during  the  heating 
mode  was  smaller  than  during  the  cooling  mode  making  the  refrigerant 
quality  variation  even  greater  in  the  indoor  coil. 

If  it  is  assumed  that  these  differences  are  due  to  deficiencies  in  the  simulation 
program  in  its  present  form,  then  it  is  concluded  that  the  condensing  heat  trans- 
fer coefficient  correlation  used  in  HPSIM  provides  a  lower  value  for  the  indoor 
coil  than  actually  exists  and  is  responsible  for  the  high  refrigerant  discharge 
pressure  obtained  in  the  computer  simulation  runs. 

6.3    VERIFICATION  OF  THE  MASS  INVENTORY  SIMULATION 

Verification  of  the  mass  inventory  method  can  be  done  on  different  levels.  One 
of  the  ways  is  to  compare  the  refrigerant  charge  in  a  system  tested  in  a  labo- 
ratory to  the  charge  computed  during  a  simulation  run.     This  was  done  for  system 
1.     The  system  simulation  resulted  in  a  charge  of  6.1  lb.  while  the  actual  heat 
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Table  4.     Cooling  Test  Simulation  Results 


SYSTEM 

TEST 

TOA 

TRA 

RHRA 

QLOAD 

ELUSE 

COP 

RMASS 

P3 

TSUP3 

P8 

°F 

°F 

Btu/h 

Watt 

Ib/h 

ps  ia 

°F 

psia 

1 

1.1 

104.9 

79.8 

.52 

33750 

5248 

1.884 

557 

86.9 

2.0 

288.0 

1 

1.2 

95.2 

79.8 

.52 

35580 

5000 

2.085 

545 

82.6 

11.0 

255.4 

1 

1.3 

82.0 

79.6 

.52 

37510 

4650 

2.364 

525 

76.6 

25.0 

214.9 

1 

1.4 

82.1 

79.5 

.19 

33690 

4572 

2.159 

506 

73.4 

1.0 

212.8 

2 

2.1 

104.9 

80.0 

,52 

36010 

5401 

1.954 

586 

90.6 

5.0 

291.8 

2 

2.2 

94.9 

80.0 

.52 

37890 

5129 

2.165 

574 

86.1 

13.8 

258.2 

2 

2.3 

82.1 

80.0 

.52 

39950 
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2.448 
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79.9 

28.0 

218.5 

2 

2.4 

82.1 

79.8 

.20 

36570 

4740 

2.261 
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77.7 

1.0 

217.4 

3 

3.1 

105.1 

79.9 

.51 

37810 

5690 

1.947 

627 

95.3 

1.0 

298.8 

3 

3.2 

95.2 

79.8 

.52 

40400 

5422 

2.184 

615 

91.1 

12.8 

265.0 

3 

3.3 

81.9 

79.9 

.51 

42940 

5067 

2.483 

598 

85.1 

2.2 

222.3 

3 

3.4 

82.2 

79.9 

.21 

38680 

4990 

2.271 

577 

81.6 

1.1 

222.0 

4 

4.1 

104.9 

80.0 

.51 

32040 

4644 

2.022 

514 

90.5 

19.8 

326.9 

4 

4.2 

95.0 

80.0 

.51 

33070 

4379 

2.213 

498 

86.9 

25.0 

291.7 

4 

4.3 

82.2 

80.0 

.51 

33970 

4009 

2.482 

466 

80.4 

33.2 

248.0 

4 

4.4 

82.1 

80.0 

.18 

32020 

3952 

2.374 

456 

77.1 

21.3 

244.0 
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Table  5.     Heating  Teat  Simulation  Results 


SYSTEM 

TEST 

TOA 

RHOA 

TRA 

QLOAD 

ELUSE 

COP 

P3 

TSUP3 

P8 

Of 

°F 

Btu/h 

Watt 

psia 

°F 

psia 

2 

2.5 

46.6 

.73 

70.2 

42170 

4776 

2.587 

70.3 

0.0 

262.9 

2 

2.6 

35.0 

.82 

69.9 

35600 

4363 

2.391 

59.8 

0.0 

239  .5 

2 

2.7 

17.1 

.70 

70.0 

24400 

3668 

1.949 

43.6 

0.0 

203.0 
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Table  6.    Discrepancy  between  Laboratory  Test  Results 
and  Computer  Simulation  Results 


TEST 

OLOAD 

ELUSE 

COP 

RMASS 

P3 

TSUP3 

P8 

% 

% 

% 

% 

ps  i 

°F 

ps  i 

1.1 

1.6 

3.0 

-1.5 

-2.6 

0.3 

0.2 

2.0 

1 .2 

-1.5 

3.0 

-4.4 

-2.6 

-0.1 

0.0 

1.1 

1.3 

-2.7 

1.9 

-4.4 

-3.8 

-1.0 

1.6 

1.0 

1 .4 

1.5 

1.4 

0.1 

-5.2 

0.4 

0.0 

2.3 

2.1 

3.4 

3.6 

0.3 

-0.5 

1.7 

2.7 

-0.4 

2.2 

0.9 

3.9 

-2.9 

-0.6 

1.1 

0.0 

1.3 

2.3 

-0.4 

3.3 

-3.5 

-0.7 

0.6 

2.9 

-1.7 

2.4 

-1.6 

3.4 

-4.8 

-1.1 

1.4 

8.4 

-0.8 

2.5 

-1.1 

0.3 

-1.4 

-1.1 

0.0 

23.0 

2.6 

0.1 

-2.4 

2.6 

-0.6 

0.0 

21.3 

2.7 

-2.2 

-7.3 

5.5 

-2.0 

0.0 

15.1 

3.1 

-0.4 

2.7 

-3.0 

-0.2 

-0.2 

0.2 

-6.4 

3.2 

-1.3 

2.8 

-4.0 

-1.1 

0.2 

0.0 

-5.3 

3.3 

0.2 

2.8 

-2.5 

-0.6 

0.0 

0.1 

-6.7 

3.4 

0.3 

2.5 

-2.2 

-2.3 

1.4 

0.0 

-2.9 

4.1 

1.4 

0.7 

0.7 

3.0 

1.2 

-0.2 

10.2 

4.2 

1.3 

0.6 

0.7 

3.1 

2.1 

0.0 

11.6 

4.3 

0.9 

0.7 

0.2 

3.8 

3.8 

0.7 

11.8 

4.4 

1.9 

-0.3 

2.2 

2.5 

1.9 

4.6 

10.4 
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pump  was  charged  with  8.3  lb  of  refrigerant.     For  the  same  system,  the  outdoor 
coil  was  weighed  during  cooling  mode  operation  giving  5.4  lb  of  refrigerant 
while  the  computed  charge  for  this  coil  was  4.8  lb.     Due  to  the  test  installa- 
tion, it  was  not  possible  to  weight  the  indoor  coil;  however,  estimates  based 
on  the  differences  in  the  above  data  indicate  that  almost  three  times  more 
refrigerant  existed  in  an  indoor  coil  (working  as  evaporator)  than  that  calcu- 
lated.    A  discrepancy  of  the  same  order  of  magnitude  for  the  Martinelle-Nelson 
method  was  reported  by  Staub  and  Zuber  [37].     Though  better  agreement  between 
a  laboratory  and  calculated  charge  would  promote  more  confidence  in  the 
Lockhard-Martinelli  method,  it  should  be  recognized  that,  as  explained  in 
section  5.6,  the  absolute  value  of  the  computed  charge  is  not  critical  since 
the  program  iterates  a  new  system  status  with  change  in  operating  conditions 
based  on  refrigerant  mass  inventory  performed  at  design  conditions.     What  is 
most  important  for  modeling  purpose  is  that  the  mass  inventory  based  on  the 
Lockhard-Martinelli  method  provides  a  means  for  calculation  of  change  of  the 
refrigerant  vapor  superheat  at  the  compressor  can  inlet  sufficiently  well  for 
accurate  heat  pump  performance  prediction. 
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7.     SUMMARY  AND  CONCLUSIONS 


The  main  goal  of  this  study  was  to  formulate  an  analytical  model  of  an  air- 
source  heat  pump  equipped  with  a  constant  flow  area  expansion  device  (capillary 
tube)  that  could  provide  accurate  heat  pump  performance  predictions. 

In  the  first  step  of  this  study,  operation  of  a  heat  pump  was  analyzed  and 
previous  work  in  steady-state  vapor  compression  cycle  modeling  was  reviewed. 
Three  mathematical  models  of  a  heat  pump  were  found  in  the  literature  and  their 
applicability  and  limitations  were  studied.     Two  of  the  models  could  only  sim- 
ulate systems  equipped  with  a  thermostatic  expansion  valve  (TXV)  which  maintains 
constant  refrigerant  vapor  superheat  at  an  evaporator  outlet.     Operation  of  a 
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heat  pump  with  a  capillary  tube  is  different  since  refrigerant  vapor  superheat 
varies  at  that  location  with  changing  operating  conditions.     The  third  heat  pump 
simulation  program,  originally  formulated  for  a  system  with  TXV,  could  accommo- 
date a  capillary  tube;  however,  without  the  ability  to  calculate  refrigerant 
vapor  superheat  (assuming  its  value  is  always  known).     Also  it  was  found,  that 
the  capillary  tube  model  utilized  in  this  program  was  adequate  only  when  a  sub- 
cooled  liquid  was  entering  the  capillary  tube  entrance  and  would  not  simulate 
the  case  of  two-phase  entering  the  tube  which  happens  particularly  during  the 
heating  mode.     None  of  the  existing  models  used  operating  conditions  as  input. 

Thus,  two  major  conclusions  were  drawn.     First  that  a  new  capillary  tube  model 
had  to  be  developed  which  would  not  carry  any  modeling  restrictions.  Second, 
some  method  should  be  formulated  for  calculation  of  refrigerant  vapor  superheat 
at  the  compressor  suction  side.     It  was  also  concluded  that  a  program  performing 
heat  pump  simulation  at  imposed  operating  conditions  would  be  desirable. 

In  the-  second  step  of  this  study,  new  logic  for  a  heat  pump  simulation  program 
was  proposed  which  would  satisfy  the  identified  needs.     The  new  logic  was  based 
on  analysis  of  heat  pump  behavior  at  changing  operating  conditions  and  allowed 
for  calculation  of  refrigerant  vapor  superheat  at  the  compressor  can  inlet  by 
the  refrigerant  mass  inventory  simulation  method.     Analysis  of  the  main  heat 
pump  components  was  performed  and  phenomena  taking  place  were  described  analy- 
tically by  fundamental  equations  or  applicable  correlations.     A  comprehensive 
literature  search  was  done  and  the  best  correlations^/  available  were  selected 
for  the  model.     These  equations  were  used  to  formulate  individual  models  of 
particular  heat  pump  components,  including  a  constant  flow  area  expansion  device 
not  available  in  the  literature.     Finally,  individual  models  were  incorporated 
into  an  overall  heat  pump  model,  HPSIM. 

Verification  of  the  model  was  conducted  against  laboratory  test  data  of  three 
heat  pump  systems  in  the  cooling  mode  and  one  system  in  both  heating  and  cooling 
mode  covering  the  whole  range  of  heat  pump  operating  conditions.     The  maximum 
difference  between  laboratory  results  and  model  predictions  was  found  to  be  3.4 
percent  for  capacity  and  5.5  percent  for  coefficient  of  performance. 

The  model  HPSIM  developed  in  this  study  is  the  only  model  in  the  open 
literature  that  is  capable  of  simulating  a  heat  pump  with  a  constant  flow  area 
expansion  device  at  imposed  operating  conditions  without  restrictions  on  refrig- 
erant state  at  any  system  location.     (A  possible  exception  is  the  assumption 
that  refrigerant  saturated  vapor  is  leaving  an  accumulator  when  liquid  is  col- 
lected on  its  bottom.     However,  this  assumption  is  felt  to  be  quite  reasonable 
and  does  not  cause  modeling  difficulties.)     The  model  having  a  demonstrated  accu- 
racy of  performance  prediction,  allows  for  a  wide  variety  of  model  applications. 
As  shown  in  the  verification  process,  the  model  was  used  to  evaluate  the  per- 
formance of  systems  with  mismatched  components.     Also  it  may  be  used  to  test 
how  certain  modifications  of  components  (e.g.,  coil  tube  diameter,  fin  thick- 
ness, fin  spacing,  coil  face  area,  capillary  tube  diameter  or  length,  etc.) 


5/  Those  which  appeared  to  have  experimental  verification  in  the  conditions 
most  similar  to  those  found  during  operation  of  a  heat  pump. 
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affect  overall  performance  of  a  heat  pump  filled  with  different  refrigerants. 
This  latter  application  has  led  to  advanced  work  which  will  require  modifica- 
tion to  the  refrigerant  thermodynamic  and  transport  property  equations  to 
include  mixtures. 

The  program  HPSIM  occupies  31000  decimal  words  of  core  memory  of  a  UNIVAC  1108 
computer.     Computing  time  depends  on  size  of  coils  and  refrigerant  state  input 
guesses.     On  average,  computing  time  should  not  exceed  200  seconds.     It  is  felt 
that  the  program  can  be  helpful  in  testing,  design,  and  development  of  heat 
pump  systems  by  providing  relief  from  some  expensive  laboratory  tests  and  by 
providing  consistently  reliable  results. 
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APPENDIX  A.     CALCULATION  OF  THE  THERMODYNAMIC  AND  TRANSPORT  PROPERTIES  OF 
REFRIGERANTS  R-12  and  R-22 


Thermodynamic  property  equations  for  refrigerants  have  been  presented  in 
numerous  publications.     The  equation  of  state  was  developed  by  Martin  and  Hou 
[43].    A  review  of  correlations  and  equations  has  been  given  by  Martin  [44], 
Constants  for  these  relations  were  published  in  [45]  and  [46]  for  R-12  and 
R-22,  ^respectively .     Downing  [47]  assembled  refrigerant  thermodynamic  property 
equations  along  with  constants  for  R-12,  R-22,  and  ten  other  refrigerants. 

Equations  presented  by  Martin  [2]  have  been  supported  by  Du  Pont  and  have 
become  the  accepted  thermodynamic  properties  of  refrigerants.     They  have  been 
used  by  other  authors  in  references  [1],  [2],  [3],   [48]  and  by  several  heat 
pump  manufacturers. 

The  four  basic  refrigerant  property  equations  used  in  the  calculations  are: 

1.  Vapor  pressure  at  saturation 

2.  liquid  density  at  saturation 

3.  pressure  -  volume  -  temperature  behavior 

4.  specific  heat  of  vapor 

Refrigerant  constants  for  these  equations  can  be  found  in  table  Al.  Other 
symbols  are  explained  with  each  equation. 

The  saturation  vapor  pressure  as  a  function  of  temperature  is : 

ln(P)  =  AG  +  ^+  CG  •  ln(T)  .+  DG  •  T  +  EG(FG  -  T)  in(FG  -  T)  (Al) 
T  T 


where      P  =  vapor  saturation  pressure  (psia) 
T  =  vapor  saturation  temperature  (R) 

The  saturated  liquid  density  equation  is: 


T,        BPL  rr,         CPL  T,         DPL  ,  T     .^PL  ,  X 

Pt   =  AL  +  BL(1  --^)         +CL(1-^)         +DL(1--^)         +EL(1--!-)  (A2) 
^  TC  TC  TC  TC 

o 

where    p-^  =  saturated  liquid  dens  ity  (Ib/ft-^) 

T    =  temperature  of  saturated  liquid  (R) 

The  constant  volume  specific  heat  at  zero  pressure  equation  is: 

Cv  =  AC  +  BC  •  T  +  CC  •  t2  +  DC  •  T^  +  J|C_  (A3) 

t2 
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where    Cv  =  constant  volume  specific  heat  at  zero  pressure  (Btu/lb  •  F) 
T  =  gas  temperature  (R) 


The  equation  of  state  is: 

p  =  Al  *  T    +  A2  +  B2  «  T  +  C2  «  exp(-AK  »  T  »  TC) 
(v  -Bl)  (v  -  Bl)2 

+  A3  +  B3  »  T  +  C3  »  exp(-AK  «  T/TC)  +  A4  +  B4  «  T  +  C4  »  exp(-AK  »  T/TC) 
(v  -  Bl)3  (v  -  Bl)^ 

+  A5  +  B5  «  T  +  C5  »  exp(-AK  »  T/TC) 
(v  -  Bl)5 

+  [C6  +  B6  «  T  +  C6  *  exp(-Al  ♦  T/TC) ] exp( -2 ALPHA  *  v)  (^4) 
[CI  +  exp(-ALPHA)] 


where      P  =  vapor  pressure  (psia) 
T  =  vapor  temperature  (R) 
V  =  vapor  specific  volume  (ft-^/lb) 
TC  =  refrigerant  critical  temperature  (R) 


Using  these  four  equations  and  other  exact  thermodynamic  relations,  a  complete 
set  of  vapor  refrigerant  thermodynamic  properties  can  be  determined  under  any 
conditions . 


The  exact  thermodynamic  relations  required  are: 

2 

Cv  =  Cv  +  /  T  •  (-^)     dv  (A5) 
aT^  V 

(IP)  2 

.    Cp  =  Cv  -  T   (A6) 

(IP) 


3v 


T 


s  =  /        dT  +  /  (||)     dv  +  Y  (A7) 

* 

i  =  /Cv'dT  +  P«v  +  /  [T(-2%     -  P]dv  +  X  (A8) 

3T  y 


where  the  new  symbols  mean 


Cv  =  vapor  specific  heat  at  constant  volume 
Cp  =  vapor  specific  heat  at  constant  pressure 

i  =  vapor  enthalpy 

s  =  vapor  entropy 
X,Y  =  constant  (see  table  Al) 
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The  derivation  of  equations  (A5),  (A6),  (A7),  and  (A8)  may  be  found  in 
reference  [49]. 

Besides  the  above  equations,  algorithms  have  to  be  used  to  calculate  the 
specific  heat  at  constant  pressure  of  liquid,  as  well  as  the  viscosity  and 
conductivity  of  the  liquid  and  vapor.     There  are  many  correlations  available 
(see  reference  [1],   [2],   [3],   [50],   [51])  to  calculate  these  refrigerant  pro- 
perties.    It  was  decided  to  use  those  proposed  by  Chi  [3],  developed  in  the 
form  of  the  5th  power  polynomial,  because  they  are  accurate  within  +0.5 
percent  and  easy  to  handle  in  one  subroutine: 


6 

PROP  =    E    A(I,J)   .  T^"l  (A9) 
3=1 


where    A(I,J)  =  twelve  coefficients  for  calculated  property  (six  for  T  <  100°F, 
J  =  1  to  6  and  six  for  T  >  100°F,  J  =  7  to  12) 
1=1  for  liquid  viscosity  (Ib/h.ft) 
1=2  for  vapor  viscosity  (lb/(h«ft)) 
1=3  for  liquid  conductivity  (Btu/(h»f t ^F) ) 
1=4  for  vapor  conductivity  (Btu/(h»f t»F)) 

1=5  for  liquid  specific  heat  at  constant  pressure  (Btu/(lb»F)) 
T  =  saturated  temperature  (F) 
PROP  =  calculated  property  at  saturated  temperature 


The  coefficients  to  be  used  in  equation  (A9)  are  given  in  table  A2.  Equation 
(A9)  allows  for  the  computation  of  the  refrigerant  property  values  at  satura- 
tion temperature.     For  those  in  the  superheated  state,  there  is  a  need  to  account 
for  the  effect  of  temperature  on  viscosity  and  conductivity.     This  can  be  done 
by  using  the  power  law,  since  viscosity  and  conductivity  may  be  said  to  be  pro- 
portional to  T*^.     The  exponent,  n,  varies  with  temperature  and  is  assumed  in 
this  document  to  equal  0.5  as  called  for  by  the  kinetic  theory  of  gases  [52], 

In  conclusion,  eight  functions  and  subroutines  are  used  in  the  HPSIM  heat  pump 
simulation  program  to  determine  the  refrigerant  thermodynamic  and  transport 
properties.     The  equations  described  in  this  appendix  are  applied  with  constants 
presented  in  tables  Al  and  A2 .     The  listing  of  these  subprograms  are  given  along 
with  other  subprograms,  in  appendix  J.     Input  and  output  data  are  explained  in 
the  form  of  comment  statements  in  each  routine.     For  a  quick  reference,  a  list 
and  application  of  refrigerant  property  routines  is  also  given  in  table  A3. 
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Table  Al.    Coefficients  and  Physical  Constants  for  Freon-12  and  Freon-22 
to  be  used  in  Equations  (Al)  through  (A8) 


Coefficients  and 

Physical  Constants  FREON-12  FREON-22 

TC  6.9330000E+02  6 .6450000E+02 

PC  5.9690000E+02  7 .2191000E+02 

VC  2.8700000E+02  3.0525000E-02 

AG  9.1835883E+01  6.7598246E+01 

BG  -7.9131381E+03  -8 .8538843E+03 

GG  -1.2471522E+01  -7 .8610310E+00 

DG  1.0892245E-02  5 .0448235E-03 

EG  0.  4.4574600E-01 

FG  0.  6.8610000E+02 

AL  3.4840000E+01  3 .2760000E+01 

BL  1.8601368E+01  5 .4634409E+01 

CL  2.1983963E+01  3 .6748920E+01 

DL  5.3341175E+01  -2 .2292560E+01 

EL  -3.1509939E+00  2 .0473289E+01 

BPL  1.  3.3333333E-01 

CPL  5.0000000E-01  6 .6666667E-01 

DPL  3.3333333E-01  1. 

EPL  2.  1.3333333E+00 

Al  8.8734000E-02  1.2409800E-01 

Bl  6.5093890E-03  2 .OOOOOOOE-03 

CI  0 .  0 . 

A2  -3.4097270E+00  -4 .3535470E+00 

B2  1.5943480E-03  2 .4072520E-02 
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Table  Al .  (Continued) 


Coefficients  and 

Physical  Constants  FREON-12  FREON-22 

C2  -5.6762767E+01  -4 .4066868E+01 

A3  6.0239450E-02  -1 .7464000E-02 

B3  -1.8796180E-05  7 .6278900E-05 

C3  1.3113990E+00  1 .4837630E+00 

A4  -5.4873700E-04  2 . 3101420E-03 

B4  0.  -3.6057230E-06 

04  0 .  0 . 

A5  0.  -3.7240440E-05 

B5  3.468834OE-09  5 .3554650E-08 

05  -2.5439070E-05  -1 .8450510E-04 
A6  0.  1.3633870E+08 
B6  0.  -1.6726130E+05 

06  0.  0. 

ALPHA  0.  5.4820000E+02 

AK  5.4750000E+00  4 .2000000E+00 

AC  8.0945000E-03  2 .8128360E-02 

BO  3.3266200E-04  2 .2554080E-04 

CO  -2.4138960E-07  -6 .5096070E-08 

DO  6.7236300E-11  0. 

EC  0.  0. 

FC  0.  2.5634100E+02 

X  3.9556551E+01  6.2400900E+01 

Y  -1.6537940E-02  -4 .5333500E-02 
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Table  A2.    Coefficients  for  Freon-12  and  Freon-22  to  be  Used  in 
Equation  (A9) 


Coefficients 


A( 

:i 

,1) 

M 

,1 

.2) 

A( 

;i 

,3) 

A( 

,1 

,4) 

Al 

;i 

.5) 

A( 

,1 

.6) 

A( 

,1 

.7) 

A( 

^1 

.8) 

A( 

,1 

.9) 

A( 

^1 

,10) 

Al 

^1 

.11) 

A( 

,1 

,12) 

A( 

2 

,1) 

A( 

.2 

,2) 

A( 

2 

,3) 

A( 

^2 

A) 

A( 

2 

,5) 

A( 

2 

►  6) 

A( 

2 

,7) 

A( 

2 

,8) 

A( 

2 

,9) 

A( 

2 

,10) 

A( 

2 

,11) 

A( 

2 

,12) 

A( 

3 

1) 

A( 

3 

.2) 

A( 

3 

.3) 

A( 

.3 

A) 

A( 

3 

,5) 

A( 

3 

,6) 

A( 

3 

,7) 

A( 

3 

,8) 

A( 

3 

,9) 

A( 

.3 

,10) 

A( 

.3 

.11) 

A( 

.3 

,12) 

A( 

.4 

.1) 

A( 

[4 

.2) 

A( 

4 

,3) 

A( 

.4 

,4) 

A( 

4 

,5) 

A( 

:4 

,6) 

A( 

> 

,7) 

Freon-12 

0.75800E+00 
0.44230E-02 
0,22659E-04 

-0.80937E-07 
0.48640E-09 

-0.38993E-11 
0.16672E+01 
0.37691E-01 
0.48998E-03 

-0.32598E-05 
0.10701E-07 
0.14069E-10 

0.26200E-01 

0.58000E-04 

0.0 

0.0 

0.0 

0.0 

0.10065E+01 
0.37971E-01 
0.53504E-03 
0.36443E-05 
0.12032E-07 
0.15478E-10 

0.49000E-01 
0.11950E-03 
0.36320E-07 
0.52080E-09 
0.31690E-11 
0.31689E-13 
0.33483E+00 
0.10635E-01 
0.14902E-03 
0.10214E-05 
0.34038E-08 
0.44300E-11 

0.43000E-02 

0.17000E-04 

0.0 

0.0 

0.0 

0.0 

0.12814E+00 


Freon-22 

0.64600E+00 
0.29194E-02 
0.12164E-04 

-0.74985E-07 
0.83951E-09 

-0.37512E-11 
0.69684E+01 
0.24319E+00 
0.35924E-02 

-0.26187E-04 
0.93884E-07 
0.13301E-09 

0.26600E-01 
0.63804E-04 
0.10761E-06 
0.32061E-08 
0.43463E-10 
0.13175E-12 
0.66330E+00 
0.25757E-01 
0.37913E-03 
0.27734E-03 
0.10053E-07 
0.14474E-10 

0.63000E-01 
0.15820E-03 
0.12289E-06 
0.17453E-08 
0.52685E-11 
0.10856E-13 
0.46705E+00 
0.19421E-01 
0.28507E-03 
0.20429E-05 
0.71963E-08 
0.99460E-11 

0.48000E-02 
0.19881E-04 
0.24815E-08 
0.28518E-09 
0.62001E-11 
0.31001E-13 
0.51539E+00 
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Table  A2 .  (Continued) 


Coefficients  Freon-12  Freon-22 

A(4,8)  0.48691E-02  0.18936E-05 

A(4,9)  0.68452E-04  0.26762E-03 

A(4,10)  0.46711E-06  0.18936E-05 

A(4,ll)  0.15484E-08  0.65891E-08 

A(4,12)  0.20038E-11  0.90041E-11 

A(5,l)  0.21700E+00  0.2710UE+00 

A(5,2)  0.14160E-03  0.24054E-03 

A(5,3)  0.64705E-06  0.38936E-07 

A(5,4)  0.55390E-08  0.23481E-07 

A(5,5)  0.13779E-10  0.97345E-10 

A(5,6)  0.17912E-12  0.44953E-12 

A(5,7)  0.15409E+01  0.49002E+00 

A(5,8)  0.64802E-01  0.83123E-02 

A(5,9)  0.91255E-03  0.13105E-03 

A(5,10)  0.62141E-05  0.96844E-06 

A(5,ll)  0.20422E-07  0.36462E-08 

A(5,12)  0.26255E-10  0.52089E-11 
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Table  A3.     Refrigerant  Property  Functions  and  Subroutines 


Property  Function 
or  Subroutine 

CPCV 


ITRPR 


SATP 


SATPR 


SATT 


SATVF 


VPSV 


VPVHS 


Application 

Calculate  specific  heat  at  constant  volume  and  at 
constant  pressure,  specific  heat  ratio,  and  sonic 
velocity  of  refrigerant  vapor  from  temperature 
and  pressure. 

Iterate  refrigerant  vapor  thermodynamic  properties 
from  given  pressure  and  specific  volume  or  enthalpy 
or  entropy. 

Calculate  refrigerant  saturation  pressure  from  given 
temperature. 

Calculate  dynamic  viscosity,  thermal  conductivity  of 
liquid  and  vapor,  and  refrigerant  liquid  specific 
heat  at  saturation  from  given  temperature. 

Calculate  refrigerant  saturation  temperature  from 
given  pressure. 

Calculate  specific  volume  of  saturated  liquid 
refrigerant  from  given  temperature. 

Calculate  specific  volume  of  vapor  from  given 
temperature  and  pressure. 

Calculate  refrigerant  vapor  thermodynamic  properties 
from  given  temperature  and  pressure. 
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APPENDIX  B.     CALCULATION  OF  PROPERTIES  OF  MOIST  AIR 


For  the  selection  of  moist  air  property  equations,  it  was  assumed  that  air 
dry  bulb  temperature,  T,  relative  humidity,  and  pressure  would  be  known  and 
the  properties  that  would  be  the  output  were:     humidity  ratio,  specific  heat 
at  constant  pressure,  gas  constant,  dynamic  viscosity,  and  thermal  conduc- 
tivity.    Presented  here  are  the  psychometric  equations  in  their  fundamental 
form.     They  are  derived  assuming  moist  air  to  be  a  mixture  of  two  independent 
perfect  gases,  so  that  the  perfect  gas  equation  of  state  and  Dalton's  rule 
can  be  applied.     The  transport  properties,  dynamic  viscosity,  and  thermo- 
dynamic conductivity  of  air  are  assumed  to  be  negligibly  affected  by  the 
moisture  content.     Correlations  for  these  two  properties  were  based  on  dry 
air  data. 

Relations  of  this  appendix  were  applied  in  air  properties  subroutine  AIRPR. 

The  humidity  ratio  is  determined  by: 

w  =  0.622  PSAT/(PATM  -  (|)  •  PSAT)  (Bl) 

where    PATH  =  atmospheric  pressure 

PSAT  =  saturation  pressure  of  water  vapor  at  temperature  T 

w  =  humidity  ratio  (lb  water/lb  dry  air) 

<j)  =  relative  humidity  (-) 

The  saturated  water  vapor  pressure,  PSAT  (psi),  is  calculated  by  the 
polynominal  approximation: 

'  For  T  >  32°F 
PSAT  =  exp(19.504  -  10.431z  -  0.2755z2  +  0.039402^) 

For  32  <  to  <  180°F 

PSAT  =  exp(13.4353  -  5.0988z  -  1.6896z2  +  0.17829z3) 

For  T  >^  180°F 

PSAT  =  exp(16.8255  -  14.213z  +  7 .5568z2  -  4.01506z3  +  0.17692z^) 
1000 

where    z  =   (B2) 

460  +  T 

The  specific  heat  at  constant  pressure  is: 

Cp  =  (Cpdry  +  0.444)7(1  +  w)  (B3) 

where    Cp        =  specific  heat  at  constant  pressure  of  moist  air  (Btu/(lb  moist 
air  •  F)) 

Cpdry  =  specific  heat  at  constant  pressure  of  dry  air  (Btu/(lb  •  F)) 
w  =  humidity  ratio  (lb  water/ lb  dry  air) 
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The  specific  heat  at  constant  pressure  of  dry  air  is  approximated  by  the 
following  polynominal: 

Cp,      =  0.2478786  -  0.4204563  •  10"^  •  TR  +  0.567857  •  10"^  •  TR^ 

^  -  0.14936056  Tr3  (B4) 

where    TR  =  T  +  460*'F  is  the  dry  bulb  air  temperature  on  the  Rankine  scale 

The  gas  constant  is: 

R  =  (53.34  +  85.76  •  w)/(l  +  w)  (B5) 

Ibf  •  ft 

where    R  =  moist  air  gas  constant  (  ) 

lb  •  R 

w  =  humidity  ratio  (lb  water/lb  dry  air) 

The  dynamic  viscosity  and  thermal  conductivity  values  are  obtained  from  the 
respective  equations: 

y  =  5.5029  •  lO'^  +  8.7157  •  lO'^  TR-2.9464  •  10-8'rR2  +  6.25  •  10~'^^TR^  (B6) 

k  =  -2.853  •  10-^  +  3.268  •  lO'^  tr  -8.253  •  10-9tr2  +  1.239  •  IQ-'^^TK^  (B7) 

where      k  =  moist  air  thermal  conductivity  (Btu/(h»F» f t )) 

TR  =  460  +  T  =  dry  bulb  air  temperature  on  the  Rankine  scale  (R) 


106 


APPENDIX  C.     CALCULATION  OF  WATER  AND  FROST  PROPERTIES 


The  equations  presented  below  are  for  calculation  of  the  properties  of  water 
and  frost  deposited  on  a  heat  pump  evaporator  outer  surface.     These  properties, 
for  specific  conditions,  are  either  solely  a  function  of  temperature  or  can  be 
assumed  to  be  constant.     The  equations  below  are  used  in  the  water  and  frost 
properties  subroutine  WATPR. 

The  following  fourth  degree  polynomial  expression  is  used  for  calculating  water 
density,  conductivity,  dynamic  viscosity  and  latent  heat  of  condensation: 

5 

PROP  =    Z     A(I)   •  T^-l  (CI) 
1=1 

where    A(I)  =  five  constants  per  calculated  property 
PROP  =  calculated  property 
T  =  water  temperature  (F) 

Constants  A(I)  for  each  property  are  given  in  table  CI.  The  specific  heat  of 
water  at  constant  pressure  is  assumed  to  be  independent  of  temperature  and  to 
have  a  constant  valu  of  Cp  =  l_Btu —  . 

lb  •  F 

The  density  of  frost,  as  proposed  in  [53],  is  calculated  by  the  following 
equation: 

Pf  =  exp(bl  +  b2  •  (TW  -  TP))  (C2) 

where    bl  =  11.9521  +  0.02422  TPR  +  35.5498  WA  -9.1742  •  10-7  VA  + 
3.1138  •  10-9  VA  •  TPR  -0.03838 
b2  =  13.1606  -  0.02133  TPR  -81.955  WA/32.018  -  TP 
TP  =  tube  temperature  (F) 
TPR  =  TP  +  460.  tube  temperature  (R) 
TW  =  water  (frost)  temperature  (F) 
WA  =  air  humidity  ratio 
VA  =  air  velocity  (ft/sec) 
P£  =  density  of  frost  (lb/ft  ) 

The  frost  conductivity,  as  proposed  in  [54],  is  calculated  by  the  equation: 

=  0.012138  +  3.8909  10"^  •  pj  +  5.1409  10"^  •  p^^  (C3) 


where    kf  =  frost  conductivity  (Btu/(h*F 'f t)) 
pf  =  frost  density  (lb/ft"') 


Frost  heat  of  sublimation,  hsugi,,  and  frost  specific  heat,  Cf ,  are  assumed  to 
have  the  following  constant  values: 

hsuBL  =  1219.0  Btu/lb 
Cf  =  0.46  Btu/(lb»F) 


107 


Table  CI.    Water  Property  Evaluation  Constants  which 
are  Used  in  Equation  (CI) 


L.  ax  cuJLa  u  ea 
Property 

A(I) 

Density 
Ib/ft^ 

Dynamic 
viscosity 
lb/(h'ft) 

Thermal 
conductivity 
Btu/(h«F'ft) 

Latent  heat 
of  condensation 
Btu/lb 

A(l) 

0.11647  E+03 

0.79422  E+03 

-0.27694 

0.31514  E+04 

A(2) 

-0.40054 

0.47589  E+01 

0.45215  E-03 

0.13714  E+02 

A(3) 

0.10815  E-02 

0.10622  E-01 

0.49008  E-05 

0.35945  E-01 

A(4) 

0.12387  E-05 

0.10416  E-04 

0.88613  E-08 

0.43525  E-04 

A(5) 

0.49002  E-09 

0.37690  E-08 

0.41387  E-11 

0.19695  E-07 
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APPENDIX  D.     COMPRESSOR  SIMULATION  SUBROUTINE  PERCOM 


Reference  Is  made  to  figure  4  where  the  configuration  of  the  main  heat  pump 
components  is  shown.     It  was  decided  to  cover  in  one  subroutine  called  PERCOM 
all  refrigerant  processes  from  point  1  to  point  10.     The  refrigerant  path 
covered  by  this  subroutine  consists  of: 

-  flow  through  a  suction  tube  1-2  and  discharge  tube  9-10 

-  flow  through  a  reversing  valve  on  suction  side  2-3  and  discharge  side 

8-9 

-  flow  through  a  compressor  3-8 

Input  data  to  PERCOM  is  detailed  in  comment  statements  inserted  in  the  program. 
Basically  it  consists  of: 

--connecting  tube  parameters 

-  four-way  valve  parameters 

-  compressor  parameters 

-  ambient  air  conditions 

-  refrigerant  state  at  the  compressor  can  inlet  (point  3) 

-  refrigerant  pressure  at  discharge  (point  6) 

Compressor  and  four-way  valve  parameters  have  to  be  found  from  a  laboratory 
test.     For  a  four-way  valve,  these  parameters  are  heat  transfer  and  pressure 
drop  parameters  and  can  be  calculated  with  aid  of  the  program  VALPAR.     In  the 
case  of  a  compressor,  parameters  to  be  found  are: 

Tie  ~  electric  motor  efficiency  versus  load 
RPM  =  electric  motor  speed  (RPM)  versus  load 

Tija  =  compressor  mechanical  efficiency 

Cg  =  compressor  effective  clearance 

rip  =  compressor  polytropic  efficiency 
CQ4,5,  CQ5J,  CQ7^8»  ^Q3,C'  CQc,A  =  ^^^^  transfer  parameters 
^^3,4 >  ^^4,5 >  ^^6,7 »  ^^7,8  ~  pressure  drop  parameters 

One  of  the  parameters,  the  compressor  mechanical  efficiency,  r\^,  is  assumed 
here  to  be  constant  and  equal  to  0.96.     The  rest  of  the  parameters  can  be  cal- 
culated by  means  of  subroutine  COMPAR.     Efficiency,  ng,  and  RPM  as  functions 
of  mechanical  load  are  both  part  of  typical  electric  motor  characteristics  for 
its  class,  (see  figures  10  and  12).     Coefficients  describing  these  curves  are 
required  as  an  input  to  COMPAR.     In  the  course  of  these  calculations,  new 
coefficients  are  computed  which  retain  the  shape  of  the  characteristics  curves, 
although  a  change  in  absolute  values  of  these  coefficients  does  occur  according 
to  test  data  supplied.     The  rest  of  the  parameters  (compressor  effective  clear- 
ance, compressor  polytropic  efficiency,  five  heat  transfer  parameters,  and  four 
pressure  drop  parameters)  are  calculated  by  COMPAR  using  test  data  and  relations 
presented  in  section  5.2.2.     Subroutine  COMPAR  is  called  by  the  main  program 
MAIN.     Logic  of  PERCOM  may  be  followed  by  looking  at  the  program  listing  with 
Inserted  comment  statements.     Based  on  a  known  refrigerant  state  at  the  can 
inlet  (point  3),  refrigerant  parameters  at  the  suction  valve  (point  5)  are 
guessed  and  since  compression  pressure  is  given,  the  compression  can  be 
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computed.     Using  the  refrigerant  mass  flow  rate  from  this  computation  and 
applying  equations  given  in  section  5.2.2,  enthalpy  and  pressure  balances  are 
conducted.     If  balance  is  not  reached,  the  guess  of  the  refrigerant  state  at 
the  suction  valve,  point  5,  is  adjusted.     Calculations  are  repeated  until 
balances  are  satisfied. 

Once  refrigerant  states  at  compressor  stations  3  to  8  and  the  mass  flow  rate 
are  known,  other  refrigerant  paths  are  computed.     Simulation  of  a  four-way 
valve  and  pipes  are  done  by  routines  VALVE4  and  PIPEPI  respectively. 

Output  from  subroutine  PERCOM  is  detailed  in  the  program  comment  statements 
and  consists  of  refrigerant  thermodynamic  properties  for  points  1  to  10,  refri- 
gerant mass  flow  rate,  compressor  RPM,  energy  consumption  rate,  electric  motor 
efficiency,  compression  efficiency,  and  volumetric  efficiency. 
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APPENDIX  E.     CALCULATION  OF  CRITICAL  PRESSURE  FOR  TWO-PHASE  FANNO  FLOW 


In  this  appendix,  the  procedure  is  explained  for  calculation  of  critical 
pressure  of  a  flow  if  mass  flow  rate  and  stagnation  enthalpy  are  given.  The 
procedure  focuses  on  the  fact  that  Fanno  flow  assumes  maximum  entropy  at  the 
flow  critical  pressure.     Hence,  the  point  of  maximum  entropy  on  the  Fanno  line 
is  being  sought  and  once  determined,  the  choking  pressure  is  found. 

The  entropy  of  two-phase  flow  is  calculated  by  the  equation: 

s  =  SL  +  x(sv  -  sl)  (El) 

where    s  =  entropy 
x  =  quality 

subscripts  L  and  V  refer  to  liquid  and  vapor  phase,  respectively. 

Since  the  entropy  of  saturated  liquid  and  of  saturated  vapor,  sl  and  sy,  are 
known  for  a  given  pressure,  the  entropy  of  two-phase  flow  can  be  found  if  the 
quality  is  known. 

The  quality  for  Fanno  flow  can  be  found  using  the  energy  equation: 

°  2 

where      G  =  mass  flux 
i  =  enthalpy 
Iq  =  stagnation  enthalpy 
V  =  specific  volume 

Two-phase  specific  enthalpy  and  specific  volume  are: 

i  =  iL  +  x(iv  -  11)  (E3) 
V  =  Vl  +  x(vv  -  Vl)  (E^) 

Substituting  and  rearranging,  the  following  quadratic  equation  can  be  obtained: 
x2+x»b+c=0  (E5) 
2(v    -v)«v     'G^+i  -i 


where    b  = 


2a  -  i  )  - 
e  =   L  o_  +  ^  2 


g2  ^    (v^  -  Vl)2 

Equations  (E5)  allows  calculation  on  the  quality  for  given  Fanno  flow  (mass  flux 
and  stagnation  enthalpy  known)  at  a  particular  pressure.     This  quality  can  then 


111 


be  used  to  find  entropy  by  equation  (El).     Thus  the  maximum  entropy  and  the 
critical  pressure  of  the  flow  can  be  found  numerically^ 

Search  for  the  maximum  entropy  along  the  Fanno  flow  path  was  tried  by  two 
different  iteration  methods.     The  Newton-Raphson  method  appeared  to  be  unstable 
for  single-precision  calculations.     Stability  of  this  method  could  be  assured 
by  performing  calculations  in  the  double  precision  mode.     Such  calculations 
consume  twice  as  much  computing  time  as  single-precision  calculations  and  would 
require  keeping  some  refrigerant  property  subroutines  as  double  precision  rou- 
tines.    It  was  felt  to  be  more  reasonable  to  use  the  binary  search  method.  The 
binary  search  method  is  not  as  efficient  as  the  Newton-Raphson  method;  however, 
it  provides  a  stable  solution  in  single-precision  calculations  for  maximum 
entropy. 
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APPENDIX  F.     CAPILLARY  TUBE  SIMULATION  SUBROUTINE  EXDEV 


The  purpose  of  the  capillary  tube  simulation  routine  is  to  calculate  the 
refrigerant  mass  flow  rate  when  capillary  dimensions,  refrigerant  inlet  state, 
and  pressure  in  the  evaporator  are  given.     Depending  on  the  above  input,  there 
are  several  possible  capillary  tube  operation  modes.     For  example,  the  refri- 
gerant can  be  either  a  single  or  a  two-phase  state  prior  to  the  tube  entrance. 
If  there  is  a  two-phase  mixture  at  the  inlet  (a  highly  undesirable  but  possi- 
ble condition  under  very  low  load  operating  conditions),  two-phase  flow  will 
exist  along  the  whole  tube  length. 

For  a  subcooled  liquid  at  the  inlet,  there  are  three  flow  alternatives: 

-  liquid  flow  only  along  the  whole  tube  length 

-  liquid  flow  in  the  front  portion  of  the  tube  and  two-phase  flow  in 
the  further  portion  of  the  tube 

-  two-phase  flow  only  along  the  whole  tube  length 

Each  of  the  above  flow  situations  can  exist  with  or  without  choking  condition 
at  the  tube  exit.     The  capillary  tube  simulation  routine  EXDEV  developed  in 
this  study  can  distinguish  among  all  of  the  above  operational  modes.  During 
the  course  of  the  calculations,  the  refrigerant  critical  pressure  is  found  by 
routine  PCHOKE  using  a  binary  search  iteration  method.     Then  equations  (38), 
(43),  and  (48)  are  used  to  compute  the  refrigerant  mass  flow  rate.  The 
density-pressure  integral  in  equation  (48)  is  evaluated  by  routine  TSIMP  based 
on  a  Simpson's  rule.     An  iteration  solution  is  required  to  obtain  the  final 
simulation  result  since  choking  pressure,  friction  factor,  fractions  of  tube 
length  with  liquid  and  two-phase  flow,  and  the  velocity  head  used  to  correct 
enthalpy  are  functions  of  refrigerant  mass  flow  rate  which  has  to  be  found. 
The  Newton-Raphson  method  is  used  in  the  iteration  for  the  final  EXDEV  solution. 
The  logic  of  the  routine  EXDEV  is  presented  in  figure  Fl.     The  input  and  output 
data  are  listed  in  comment  statements  in  this  routine. 


If3 


f  INPUT: 

EXPANSION  DEVICE  LENGTH  &  INNER  DIA. 
REFRIG.  ENTHALPY  BEFORE  TUBE  ENTRANCE 
REFRIG.  PRESSURE  BEFORE  TUBE  ENTRANCE 
REFRIG.  PRESSURE  IN  EVAPORATOR 
GUESSED  REFRIG.  MASS  FLOW  RATE 


FIND  ALL  REFRIG.  STATE  PARAMETERS 
BEFORE  EXPANSION  DEVICE  ENTRANCE 


(liquid  at  inlet) 


ASSIGN  REFRIG.  SAT.  PRESSURE 
AT  TUBE  OUTLET 

FIND  MASS  FLOW  RATE  Ml 


1  FIND  PFLASH 

PFLASH  =  PO 

r 

FIND  PEXIT  1 

FIND  PI 


NO 


CALC.  LENGTH  OF  TUBE 
WITH  SU6C00LED  LIQUID 

FIND  LENGTH  OF  TUBE 
WITH  2-PHASE  MIXTURE 


1*  PFLASH^ 


YES 
|2-phase 
flow  only 
in  the  tube] 


ADJUST  Ml 


CALC.  REFRIG.  MASS  FLOW  RATE  M2 
FOR  2-PHASE  PORTION  OF  TUBE 


Figure  Fl.     Logic  of  constant  flow  area  expansion  device  simulati 
program  EXDEV 
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APPENDIX  G.     EVAPORATOR  AND  CONDENSER  SIMULATION  SUBROUTINES  EVPHXM  AND  CNDNEW 


The  tube-by-tube  approach  to  coil  simulation  has  been  explained  in  section  5.4 
where  applicable  heat,  mass,  and  momentum  transfer  equations  have  also  been  pre- 
sented.    Though  the  evaporator  and  condenser  simulation  programs  require  some- 
what different  heat  transfer  relations,  their  logic  is  the  same  as  shown  in 
figure  Gl. 

The  iilput  to  coil  simulation  program  is  coil  design  data,  refrigerant  flow 
distribution  sequence,  and  air  and  refrigerant  inlet  conditions.     During  the 
first  bank,  of  tubes'  performance  calculations,  the  air  parameters  upstream  of 
each  tube  row  are  estimated.     The  program  then  computes  the  heat  transfer  rate 
and  pressure  drop  tube-by-tube  from  the  inlet  to  the  exit  of  the  coil.  Mean 
refrigerant  properties  in  each  tube  are  used  to  calculate  heat  transfer  rates 
and  pressure  drops  applying  equations  described  in  section  5.4.    An  iterative 
procedure  is  used  for  each  tube  since  only  refrigerant  conditions  at  tube  inlet 
are  known  when  the  calculations  start,  and  initially  the  average  refrigerant 
propetties  in  a  tube  have  to  be  assumed  to  be  those  at  the  inlet  conditions. 
Upstream  air  temperature  and  humidity  ratio  may  not  be  accurately  known  at  the 
time  tube  performance  is  being  calculated,  so  a  second  iterative  loop  is  used 
to  update  air-side  data  after  each  round.     It  should  be  realized  that,  while 
in  the  condenser  case,  update  of  air-side  data  means  update  of  air  temperature 
and  humidity  ratio;  for  the  evaporator  it  also  includes  liquid  (frost)  layer 
thickness.     When  refrigerant  enthalpy  values  at  the  coil  exit  obtained  from  two 
consecutive  loops  are  within  the  imposed  tolerance,  the  calculations  of  coil 
performance  are  completed  and  results  printed. 

Both  the  evaporator  coil  simulation  model  EVPMXM  and  condenser  coil  simulation 
model  CNDNEW  comply  with  the  logic  and  equations  previously  described^/ .  Input 
and  output  data  for  these  programs  are  detailed  in  comment  statements  at  the 
beginning  of  each  subroutine.     Other  inserted  comment  statements  along  with  the 
flow  chart  (figure  Gl)  facilitate  the  understanding  of  both  program's 
organization. 


6/ 

For  reasons  of  organization  of  the  overall  heat  pump  simulation  program, 
EVPHXM  performs  evaporator  calculations  backwards  from  the  evaporator 
exit  to  the  inlet. 
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QarT]) 


'INPUT: 

COIL  GEOMETRY 

REFRIGERANT  FLOW  DISTRIBUTION  i 
AIR  INLET  CONDITIONS 
REFRIGERANT  INLET  CONDITIONS 

I   


ESTIMATE  AIR  PARAMETERS 
UPSTREAM  EACH  TUBE  DEPTH  ROW 


PRINT  RESULTS 


Figure  Gl.    Flow  chart  for  coil  simulation  program 
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APPENDIX  H.     THE  LOGIC  FOR  THE  MAIN  PROGRAM,  MAIN 


The  principle  of  the  logic  for  the  heat  pump  model  has  been  explained  in  section 
5.1.2.     It  is  based  on  three  balances: 

-  enthalpy  balance 

-  pressure  balance 

-  mass  balance 

The  enthalpy  balance  and  pressure  (mass  flow  rate)  balance  are  interdependent 
and  have  to  be  found  in  a  simultaneous  iteration  process.     Their  solution  pro- 
vides refrigerant  properties  in  the  system  which  are  required  for  mass  inven- 
tory calculations.     The  balances  are  performed  in  the  manner  as  shown  in  figure 
HI,  which  presents  the  logic  developed  during  this  study  and  contained  in  the 
main  program  MAIN. 

The  objective  of  the  logic  is  to  iterate  refrigerant  states  at  key  locations 
of  a  heat  pump  for  given  outdoor  and  indoor  conditions.     The  addresses  of  the 
key  locations  used  in  the  program  MAIN  are  consistent  with  those  marked  in 
figures  4  and  HI: 

1  -  evaporator  exit,  suction  pipe  inlet 

2  -  suction  pipe  outlet,  low  pressure  4-way  valve  inlet 

3  -  low  pressure  4-way  valve  exit,  compressor  can  inlet 

4  -  inside  compressor  can 

5  -  compressor  cylinder  at  suction 

6  -  compressor  cylinder  at  discharge 

7  -  discharge  manifold 

8  -  compressor  can  exit,  high  pressure  4-way  valve  inlet 

9  -  high  pressure  4-way  valve  exit,  discharge  pipe  inlet 

10  -  discharge  pipe  outlet,  condenser  inlet 

11  -  condenser  outlet,  liquid  line  inlet 

12  -  liquid  outlet,  expansion  device  outlet 

13  -  evaporator  inlet 

The  thermodynamic  processes  shown  in  figure  HI  and  figure  5  are  covered  by 
five  subroutines:     PERCOM,  CNDNEW,  LQLINE,  EXDEV  and  EVPHXM.     The  iteration 
process  starts  with  a  guessed  refrigerant  pressure  and  vapor  superheat  at  the 
compressor  can  inlet  and  with  guessed  compressor  discharge  pressure  (points 
3  and  6).     From  these  data,  the  compressor  performance  simulation  subroutine 
PERCOM  computes  the  refrigerant  mass  flow  rate  through  the  compressor  and 
refrigerant  parameters  from  point  1  to  10.     The  condenser  and  liquid  line 
subroutines  CNDNEW  and  LQLINE  are  called  next  to  perform  calculations  which 
yield  refrigerant  states  at  point  11  and  12.     Subsequently,  the  expansion 
device  simulation  program  EXDEV  computes  refrigerant  mass  flow  rate  through 
the  flow  restrictor.     At  this  point,  the  mass  flow  rate  through  the  compressor, 
RMASS,  and  mass  flow  rate  through  the  expansion  device,  XMASS,  are  compared 
and  mass  flow  balance  is  sought.     This  is  done  by  resuming  calculations  by 
PERCOM,  CNDNEW,  LQLINE,  and  EXDEV,  holding  constant  refrigerant  parameters  at 
the  compressor  can  inlet  (point  3)  and  adjusting  refrigerant  saturation 
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temperature  at  discharge  (point  6)  until  an  appropriate  saturation  temperature 
at  point  6  is  found  at  which  RMASS  and  XMASS  are  equal. 

The  evaporator  performance  simulation  routine,  EVPHXM,  was  chosen  in  this  logic 
to  be  run  after  the  mass  flow  balance  between  the  compressor  and  expansion 
device  has  been  found.     (EVPHXM  requires  the  most  computing  time  of  all  the 
subroutines.)     The  subroutine  EVPHXM  computes  evaporator  performance  in  a  back- 
ward scheme  from  the  actual  flow  direction  with  the  refrigerant  state  at  point 
1  and  mass  flow  rate  as  input,  yielding  the  refrigerant  state  at  point  13  as 
output.     As  the  thermodynamic  processes  in  the  liquid  line  and  flow  restrictor 
are  adiabatic,  enthalpy  at  point  12,  H12 ,  should  be  equal  to  that  at  evaporator 
inlet,  H13.     If  they  are  not  equal,  the  refrigerant  state  at  point  3  (i.e., 
saturation  temperature)  is  adjusted  and  all  calculations  have  to  be  repeated 
from  the  beginning.     (The  other  parameter  defining  the  refrigerant  state  at 
point  3,  the  refrigerant  superheat,  is  left  unchanged.) 

Once  an  enthalpy  balance  and  mass  flow  balance  are  established,  two  out  of  three 
refrigerant  parameters  guessed  at  the  beginning  are  temporarily  obtained.  To 
verify  the  third  guessed  parameter,  i.e.,  the  refrigerant  vapor  superheat  at 
the  compressor  can  inlet,  the  refrigerant  mass  inventory  in  the  system  is  used. 
Each  time  it  is  performed,  it  provides  a  result  based  on  refrigerant  states  in 
the  system  found  after  solving  enthalpy  and  pressure  balances  with  the  assumed 
vapor  superheat  at  point  3.     If  the  amount  of  refrigerant  obtained  from  mass 
inventory  calculations  is  smaller  than  the  refrigerant  put  into  the  system,  the 
superheat  guess  has  to  be  decreased  and  all  calculations  have  to  be  repeated 
from  the  beginning.     If  the  superheat  is  equal  to  zero  and  the  mass  inventory 
result  is  still  smaller  than  refrigerant  mass  input,  it  means  that  some  of  the 
refrigerant  is  stored  in  the  form  of  liquid  in  the  accumulator  and  saturated 
vapor  is  entering  the  compressor.     This  situation  or  the  equality  between  the 
mass  inventory  results  at  design  and  actual  operating  conditions  marks  an  end 
of  the  iteration  process  and  the  results  are  printed. 

In  the  course  of  the  computing  process,  the  main  program  gathers  information 
about  the  heat  pump  components'  performance,  updated  at  each  iteration  loop, 
and  applies  them  to  anticipate  changes  of  some  parameters  caused  by  a  change 
in  some  state  property.     This  allows  for  the  iteration  process  to  converge 
faster  with  each  iteration  loop.     The  Newton-Raphson  iteration  procedure  is 
used  in  the  solution  of  the  main  program. 
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APPENDIX  I.     PROGRAM  USER'S  GUIDE 


HPSIM  is  written  in  Fortran  IV  and  makes  use  of  standard  Fortran  mathematical 
functions.     It  is  built  around  the  main  program,  MAIN,  and  41  subprograms  for 
heat  pump  component  simulation,  and  heat  transfer,  fluid  mechanics,  and  fluid 
property  calculation.     Relations  used  in  the  heat  pump  component  simulation 
subprograms,  and  fluid  mechanics  and  heat  transfer  routines  are  described  in 
the  main  text  of  this  report.     Equations  used  in  the  refrigerant,  moist  air, 
water,  and  frost  property  routines  are  described  in  appendices  A,  B,  and  C. 

All  subprograms  and  routines  used  in  heat  pump  simulation  program  are  listed  in 
alphabetic  order  in  table  II,  with  a  short  statement  of  their  application.  The 
program  HPSIM  occupies  31000  decimal  words  of  core  memory  of  Univac  1108 
computer . 

I.l     INPUT  DATA  CODING 

There  are  five  categories  of  input  values  for  the  HPSIM  runs,  namely: 

-  heat  pump  data; 

-  refrigerant  physical  constants; 

-  run  controlling  data; 

-  indoor/outdoor  air  conditions; 

-  refrigerant  parameters  (guesses). 

Refrigerant  physical  constants  and  heat  pump  design  data  are  read  into  the 
program  from  a  disk  file.     File  7  and  file  8  are  used  in  reading  statements 
which  input  these  data  for  the  refrigerant  properties  and  heat  pump  design  data, 
respectively.    The  physical  constants  for  refrigerants  12  and  22  in  proper  input 
format  are  given  in  tables  12  and  13.    The  heat  pump  design  data  input  code  is 
presented  in  table  14.     Included  are  the  Fortran  sjnnbols  with  their  explanation. 
One  example  of  heat  pump  input  data  in  original  format  is  given  in  table  15. 

Run  control  data,  indoor/outdoor  air  conditions,  design  data,  and  refrigerant 
property  guesses  are  input  at  the  terminal.     They  are  clearly  solicited  by  the 
program  and  responses  have  to  be  given  as  follows: 

1.     Request  PERFORMANCE(l)  or  COMPRESSOR  PARAMETERS( 2) ,  IPTP  = 

Response:     1  for  heat  pump  performance  evaluation 


or 


2  for  compressor  parameters  evaluation 


2. 


Request : 


ANSWER  1  FOR  YES  OR  0  FOR  NO 

DO  YOU  WANY  ANY  INPUT  DATA  PRINTED?     LPF  = 


Response : 


0  for  no  input  data  print 


or 


1  for  input  data  printout  (request  for  specification  of 
desired  data  will  follow) 
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3.  Request:       NSYS  =  1  FOR  HEATING,  NSYS  =  2  FOR  COOLING 

LONG  =  1  IF  REFRIGERANT  CHARGE  BALANCE  REQUIRED, 
OTHERWISE  PUT  LONG  =  0 

Response:     as  explained  in  the  request 

4.  Request:      TG3,  TSUP3,  TG6 

Response:     TG3  =  saturation  temperature  of  refrigerant  vapor  at 
compressor  can  inlet  (F) 

TSUP  3  =  refrigerant  vapor  superheat  at  compressor  can 
inlet  (F) 

TG6  =  saturation  temperature  of  refrigerant  vapor  in 
compressor  cylinder  after  compression  (all  guesses, 
subject  to  iteraction) 

5.  Request:       POA,  TOA,  RHOA,  PRA,  TRA,  RHRA  =  ? 
Response:     POA  =  outside  air  pressure  (psia) 

TOA  =  outside  air  temperature  (F) 

RHOA  =  outside  air  relative  humidity  (decimal  fraction) 
PRA  =  indoor  air  pressure  (psia) 
TRA  =  indoor  air  temperature  (F) 

RHRA  =  indoor  air  relative  humidity  (decimal  fraction) 

6.  Request:       I  =  999  for  DISCONTINUATION? 

Response:  999  to  terminate  HPSIM  runs,  otherwise  enter  any  other 
number,  then  system  will  repeat  requesting  data  from  1 
to  5  for  the  next  HPSIM  run. 

All  input  data  symbols  follow  general  Fortran  practice.     Variables  beginning 
with  I  to  N  are  integers,  and  all  other  letters  represent  real  numbers.  Free 
field  input  format  is  used  for  all  input  variables. 

1.2     OUTPUT  DATA  CODING 

Output  data  consists  of  refrigerant  thermodynamic  states  at  13  heat  pump  key 
locations,  identified  on  figure  4,  along  with  results  describing  the  heat  pump 
performance.     Symbols  used  for  output  data  have  the  following  meaning: 

T,P,H,S,X  -  refrigerant  temperature  (F),  pressure  (psia),  enthalpy  (Btu/lb) , 
entropy,  (Btu/(lb  •  F),  quality  (-) 
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Table  II.     Programs  Used  in  a  Heat  Pump  Simulation  Model  HPSIM 


PROGRAM 

AIRHT 

AIRPR 

CALMAS 

CNDHTC 

CNDNEW 

COMPAR 

CPCV 

DISFLW 

DYNADP 

DYNDP2 

EVPDP 

EVPHTl 

EVPHXM 

EXDEV 

FEEDP 

FINEFF 

FLASH 

ITRPR 

LQLINE 
MAIN 

OVLHTC 


PURPOSE 

Calculate  air-side  finned  tube  heat  transfer  coefficient 
Calculate  moist  air  properties 
Calculate  mass  of  refrigerant  in  a  coil 

Calculate  condensation  heat  transfer  coefficient  inside  a 
tube 

Simulate  condenser  performance 
Calculate  compressor  performance  parameters 
Calculate  specific  heat  of  vapor  refrigerant 
Determine  refrigerant  flow  distribution  in  a  coil 
Calculate  dynamic  pressure  drop  for  single-phase  flow 
Calculate  dynamic  pressure  drop  for  two-phase  flow 
Calculate  frictional  evaporation  pressure  drop 
Calculate  evaporation  heat  transfer  coefficent 
Simulate  evaporator  performance 
Simulate  expansion  device  performance 
Calculate  correction  factor  for  2-phase  pressure  drop 
Calculate  fin  efficiency 

Calculate  flashing  temperature  and  pressure  for  Fanno  flow 

Iterate  refrigerant  vapor  thermodynamic  properties  from  given 
pressure  and  other  than  temperature  property 

Calculate  pressure  drop  in  a  liquid  line 

Read  input  data,  iterate  heat  pump  simulation  solution,  print 
results 

Calculate  overall  heat  transfer  coefficient  for  dry  finned 
tube 
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PROGRAM  PURPOSE 

PCHOKE  Calculate  critical  pressure  for  two-phase  Fanno  flow 

PERCOM  Simulate  compressor  performance 

PIPEPI  Simulate  pipe  performance 

PMASS  Calculate  refrigerant  mass  in  a  tube 

RHOM  Calculate  two-phase  mixture  density  on  Fanno  flow  path 

SATP  Calculate  refrigerant  saturation  pressure 

SATPR  Calculate  refrigerant  properties  at  saturation 

SATT  Calculate  refrigerant  saturation  temperature 

SATVF  Calculate  specific  volume  of  saturated  refrigerant  liquid 

SFANNO  Calculate  two-phase  mixture  entropy  for  Fanno  flow 

SPriDP  Calculate  frictional  single-phase  pressure  drop  in  a  tube 

SPHDPl  Calculate  frictional  single-phase  pressure  drop  in  a  tube 

SPHTC  Calculate  single-phase  heat  transfer  coefficient 

TSIMP  Integrate  using  Simpson's  rule 

UPWAO  Calculate  overall  conductance  for  wet  finned  tube 

VALPAR  Calculate  reversing  valve  performance  parameters 

VALVE4  Simulate  reversing  valve  performance 

VPSV  Calculate  refrigerant  vapor  specific  volume 

VPVHS  Calculate  refrigerant  parameters  at  and  above  saturation 

WATPR  Calculate  water  and  frost  properties 

XTTCAL  Calculate  Lockhard-Martinelli  parameter 
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Table  12.     Date  File  DATA12  Containing  Physical  Constants  for  Refrigerant  12 
in  the  Input  Format  to  Program  HPSIM 


0$0$Q$*DATA12(  1  ) 


1  FREON-12 

2  6. 9330O00E+02 , 5 . 969000CE+02 , 2 . 8700000E  -  02 

3  4.5967000E+02, 1 .8505300E-01 . 2 . 7 1 828 1 8E+00 

4  9. 1335883E+01 . -  7 . 91 3 1 38 1 E+03 . -1  .2471522E+01 

5  1 .0892245E-02.0. ,0. 

6  120. .312. 

7  3.484000QE+01 , 1 .8691368E+01 .2. 1983963E+01 

8  5.3341 175E+01 . -3. 1509939E+00, 1 . 

9  5.0000000E-01 .3.3333333E-01 ,2. 

10  8.8734000E-02.6.5093890E-03 ,0. 

1 1  -3.4097270E+00, 1 . 5943480E -  03 , -5 . 6762767E+01 

12  6.0239450E-02, -1 . 87961 80E- 05. 1.31 13990E+C0 

13  -5.4873700E-04,D. .0. 

14  0.  .3.4688340E-09. -  2 . 5439070E -  05 

15  0. .0. .0. 

16  3. . 5. 4750000E+00 

17  B.0945000E-C3.3.32(p6200E-04,  -2.4138960E-07 

18  6. 7236300E- 1 1 .0. .0. 

19  3.9556551 E+01 .- 1 .6537940E-02 

20  0.75800E+00, -0.44230E-02.0.22659E-04 

21  -0.80937E-07.0.48640E-09. -0.38993E- 1 1 

22  0. 16672E+01 . -0.37691 E -01 .0.4a998E-03 

23  -0. 32598E-05.0. 10701 E -07 , -0. 14059E- 1 0 

24  0.26200E-01 .0.58000E-04.0.0 

25  0. .0. .0. 

26  -0. 10065E+01 .0. 37971 E-01 , -0.53505E-03 

27  0.56446E-05. -0. 12032E-07.0. 15478E-10 

28  0.49000E-01  . -0. 1 1 95GE -  03 . 0 . 36320E -  07 

29  0.52080E-09. -0. 31690E- 1 1 , -0 .316B9E- 1 3 

30  0.33483E+00, -0. 10635E-01 .0. 1 4902 E- 03 

31  -0. 1021 4E -05.0 . 34028E-0B . -0.44300E- 1 1 

32  0.43000E-02.O. 17000E-04.0.0 

33  0 . . 0 . , 0 . 

34  -0. 12814E+00, 0.48691 E-02 , -0.68452E-04 

35  0.4671 1E-06, -0. 1 5484E -  08 . 0 . 20038E - 1 1 

36  0.21700E+00.0. 1 41 60E -  03 , 0 . 64705E -  05 

37  0. 55390E-08. -0. 13779E- 1 0,-0.17912E-12 

38  -0. 15410E+01 .0.64802E-01 , -0 . 91 255E-03 

39  0.62141E-05, -0.20422E-G7.0.26255E-10 
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Table  13.     Data  File  DATA22  Containing  Physical  Constants  for  Refrigerant 
22  in  the  Input  Format  to  Program  HPSIM 


Q$Q$Q$*DATA22( 1 ) 


1  FREON-22 

2  6.645000CE+02 ,7.2191 OOOE  +  02 , 3 . 0525000E -  02 

3  4.5967000E+02.1  .8505300E-01  ,  2 . 7 1 828 1 8E  +  00 

4  6. 7598246E+01  , - 8 . 8538843 E+03 , -  7 . 861 031 OE+00 

5  5.0448235E-03,4.4574700E-01 , 6 . 86 1 OOOOE+02 

6  140. ,240. 

7  3. 2760000E+01 . 5 . 4634409E+0 1 . 3 . 6748920E+01 

8  -2 . 2292550E+01 , 2 . 0473289E+01 . 3 . 3333333E -01 

9  6.6656667E-01 . 1 . , 1 . 3333333E+00 

10  1 .2409800E-01  . 2  .  OOOOOOOE - 03 . 0 . 

11  -4.3535470E+00.2.4072520E-03, -  4 . 4066868E+01 

12  -1 .7464000E-02.7.6278900E-05, 1 .4837630E+00 

13  2.3101420E-03, -  3 . 6057230E -  06 , 0 . 

14  -3 . 7240440E-05 , 5 . 3S54650E -  08 , M  . 845051 0E-C4 

15  1 .3633870E+08, -1 . 6726 1 20E+05 , 0 . 

16  5.4820000E+02,4.2000000E+00 

17  2.8128360E-02.2.255403GE-04, -6 . 5096070E-08 

18  0. ,0. , 2 . E734100E+02 

19  6.2400900E-i-01  , -4.53335G0E-02 

20  0 . 64600E+00 , -  0 . 29 1 94 E -  02 , 0 . 1 2 1 64 E -  04 

21  -0 . 74985E-07 .0 . 83951 E - 09 , -  0 . 375 1  2 E -  1  1 

22  0.69684E01 , -  0 . 243 1 yEOO , 0 . 35924E -  02 

23  -0. 261 87E- 04.0. 93884 E- 07 .-0.1 3301 E -09 

24  0.26600E-01  ,  0  .  6380-JE  -  04  ,  0  .  1  0761  E- 06 

25  -0. 32061 E- 08. 0 . 43463E- 1 0.-0.13175E-12 

26  -0.66  3  30E00,0.25757E-01 . -0.37913E-03 

27  0.27734E- 05,-0. 1 0053 E- 07,0. 14474E-10 

28  0.63000E-01 , -0. 1 5820 E- 03 ,-0.1 2289E-06 

29  0. 17453E-08, -0.52635E-1 1,-0. 10857E- 13 

30  -0.46705E00,0. 19421E-01 , -0. 28507 E- 03 

31  0.204  2  9E-C5. -  0 . 7 1 963 E -  03 , 0 . 99460E - 1 1 

32  0. 48000E-02 ,0. 1 9881 E- 04 , 0. 2481 5E-08 

33  0. 2851 8E-09. -0 . 62001 E-11 .0.31001E-13 

34  0.51539E00. -0. 1 8601 E- 01 , 0. 26762 E- 03 

35  -0. 1 8936E- 05,0. 65891 E -03 . -0. 90041 E- 1 1 
35  0. 271 OOE+00,0. 24054E-03 . 0. 38936E-07 

37  0.23481E-07.-0.97345E-10.0.44953E-12 

38  0.49002E00. -  0 . 831 23E -  02 . 0 . 1 3 1 05E -  03 

39  •0.95084E-06,0.35462E-08, -0.52089E-1 1 
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Table  14.     Heat  Pump  Input  Data  Code  to  Program  HPSIM 


All  input 
same  line 

Line  1 . 


line  2; 


Line  3; 


Line  4; 


Line  5; 


Line  6 ; 


Line  7 


data  are  in  FORTRAN  free  field  input  format  with  data  values  on  the 
separated  by  commas. 

ATITLE 

ATITLE  =  title,  maximum  80  characters 
ELEFUL,  SWPVOL,  ETAPLY,  CLREFF 

ELEFUL  =  compressor  motor  energy  input  rate  at  max.  rated  load 
(kW) 

SWPVOL  =  compressor  swept  volume  per  revolution  (in^) 
ETALPY  =  compressor  polytropic  efficiency  (-) 

CLREFF  =  compressor  clearance  volume  as  fraction  of  stroke  volume 
(-) 


EMETA( I 
EMETA( I 


EMETA(I 
EMETA(I 


EMOPT(I 
EMOPT(I 
EMOPT(I 
EMPOT(I 
EMRPM( I 
EMRPM( I 


.  I  =  1,5 

=  compressor  motor  efficiency  in  fraction  at  fraction  of 
full  load  specified  by  EMPOT(I),  I  =  1,5  (-) 

I  =  6,11 

=  compressor  motor  efficiency  in  fraction  at  fraction  of 
full  load  specified  by  EMOPT(I),  I  =  6,11  (-) 

,  I  =  1,5 

=  compressor  motor  full  load  fraction  (-) 
,  I  =  6,11 

=  compressor  motor  full  load  fraction  (-) 
,  I  =  1,6 

=  coefficient  for  compressor  motor  RPM  calculations  (-) 


Line  8:       CPC34,  CPC45,  CPC67,  CPC78 

CPC34  =  pressure  drop  parameter  at  compressor  can  inlet 
((Ibf  •  h2)/(lb   •  in2  •  ft^)) 

CPC45  =  pressure  drop  parameter  at  compressor  suction  valve 
((Ibf  •  h2)/(lb   •  in2  .  ft^)) 
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CPC67  =  pressure  drop  parameter  at  compressor  delivety  valve 
((Ibf  •  h2)/(lb  •  in2  .  ft^)) 

CPC78  =  pressure  drop  parameter  at  compressure  can  exit 
((Ibf  •  h2-2)/(ib  .  In2  .  ft2-8)) 

Line  9;  '    CQC4C,  CQCCOA,  CQC45 ,  CQC67 ,  CQC78 

CQC4C  =  parameter  for  compressor  can  wall  -  refrigerant  vapor  heat 
transfer  (ft0«2) 

CQCCOA  =  parameter  for  compressor  can  -  ambient  air  heat  transfer 
(Btu/h  •  f1*333) 

CQC45  =  suction  valve  heat  transfer  parameter  (ft^*2) 

CQC67  =  delivery  valve  heat  transfer  parameter  (ft^*2) 

CQC78  =  heat  transfer  parameter  at  can  exit  (ft^*2) 

Line  10;     CQ,  CPDR 

CQ  -  parameter  for  4-way  valve  heat  transfer 

(ftO.2) 

CPDR  =  pressure  drop  parameter  for  4-way  valve 
((Ibf  •  h2)/(lb  •  in2  .  ft^)) 

Line  11;     ACCH,  ACCD,  VCAN,  REFIN 

ACCH  =  accumulator  height  (in) 

ACCD  =  accumulator  inner  diameter  (in) 

VCAN  =  volume  of  compressor  can  filled  by  refrigerant  (in^) 

REFIN  =  refrigerant  charge  (lb) 
Line  12;     NDEP(l),  NROW(l) 

NDEP(l)  =  number  of  indoor  coil  tube  depth  rows  (-) 

NROW(l)  =  rumber  of  tubes  per  indoor  coil  depth  row  (-) 
Line  13;     DI(1),  D0(1),  DT(1),  RPCH(l),  DPCH(l),  WIDTH(l) 

DI(1)  =  inner  diameter  of  indoor  coil  tubes  (in) 

D0(1)  =  outer  diameter  of  indoor  coil  tubes  (in) 

DT(1)  =■  indoor  coil  fin  tip  diameter  (in) 
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RPCH(l)  =  pitch  between  tubes  of  the  same  depth  in  indoor  coil  (in) 

DPCH(l)  =  depth  pitch  for  indoor  coil  tubes  (in) 

WIDTH(l)  =  indoor  coil  width  (equal  tube  length)  (in) 

Line  14:     FPCH(l) ,  FTK(l),  FMK(l),  TMK(l) ,  AMAS(l) 

FPCH(l)  =  indoor  coil  fin  pitch  (in) 

FTK(l)  =  indoor  coil  fin  thickness  (in) 

FMK(l)  =  indoor  coil  fin  material  thermal  conductivity 
(Btu/(ft  •  h  •  F)) 

TMK(l)  =  indoor  coil  tube  material  thermal  conductivity 
(Btu/(ft  •  h  •  F)) 

AMAS(l)  =  air  mass  flow  rate  through  indoor  coil  (Ib/h) 

Line  15;     CONST(l),  CPOW(l),  ANGLE(l) 

CONST(l)  =  constant  for  air  side  heat  transfer  correlation  for 
indoor  coil  equal  to  0.134  (-) 

GPOW(l)  =  constant  for  air  side  heat  transfer  correlation  for 
indoor  coil  equal  to  0.681  (-) 

ANGLE(l)  =  angle  between  indoor  coil  face  and  air  streamlines  (rad) 


Line 

16: 

EIDFAN 

EIDFAN  =  indoor  fan  energy  input  rate 

(kW) 

Line 

17: 

NREPTI 

NREPTI  =  number  of  repeating  sections 

in  indoor 

Line 

18: 

NTUBE(1,I)  I  =  1,5 

NTU1JE(1,1)  =  number  of  tubes  in  first  row  in  each  section  of  indoor 
coil  (-) 

NTUBE(1,2)  =  number  of  tubes  in  second  row  in  each  section  of  indoor 
coil  (-) 

NTUBE(1,3)  =  number  of  tubes  in  third  row  in  each  section  of  Indoor 
coil  (-) 

NTUBE(1,4)  =  number  of  tubes  in  fourth  row  in  each  section  of  indoor 
coil  (-) 
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NTUBE(1,5)  =  number  of  tubes  in  fifth  row  in  each  section  of  indoor 
coil  (-) 


Line  19;     IFR0M(1,I),  I  =  1,10 

IFR0M(1,1)  =  number  of  tube  of  indoor  coil  from  which  tube  1  receives 
refrigerant  when  indoor  coil  works  as  evaporator  (-) 

IFR0M(1,2)  =  number  of  tube  of  indoor  coil  from  which  tube  2  receives 
refrigerant  when  indoor  coil  works  as  evaporator  (-) 

IFR0M(1,3)  = 
IFR0M(1,9)  = 

IFR0M(1,10)  =  number  of  tube  of  indoor  coil  from  which  tube  10  receives 
refrigerant  when  indoor  coil  works  as  evaporator  (-) 

Line  20;     IFR0M(1,I),  I  =  11,20 

IFR0M(1,I)  =  number  of  tube  indoor  coil  from  which  tube  I  receives 
refrigerant  when  indoor  coil  works  as  evaporator  (-) 


Line 

21; 

IFR0M(1,I), 

I 

21,30 

Line 

22; 

IFR0M(1,I) , 

I 

31,40 

Line 

23: 

IFR0M(1,I) , 

I 

41,50 

Line 

24; 

IFR0M(1,I) , 

I 

51,60 

Line 

25; 

IFR0M(1,I), 

I 

61,70 

Line 

26; 

IFR0M(1,I) , 

I 

71,80 

Line 

27; 

IFR0M(1,I) , 

I 

81,90 

Line 

28: 

IFR0M(1,I) , 

I 

91,100 

Line 

29; 

IFR0M(1,I) , 

I 

101,110 

Line 

30; 

IFR0M(1,I) , 

I 

111,120 

Line 

31; 

IFR0M(1,I) , 

I 

121,130 

IFR0M(1,I)  =  number  of  indoor  coil  from  which  tube  I  receives 

refrigerant  when  indoor  coil  works  as  evaporator  (-) 
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Line  32;     NDEP(2),  NR0W(2) 

NDEP(2)  =  number  of  outdoor  coil  tube  row  depths  (-) 

NR0W(2)  =»  number  of  tubes  per  outdoor  coll  depth  row  (-) 
Line  33;     DI(2),  D0(2),  DT(2),  RPCH(2),  DPCH(2),  WIDTH(2) 

DI(2)  =  Inner  diameter  of  outdoor  coll  tubes  (In) 

D0(2)  =«  outer  diameter  of  outdoor  coll  tubes  (In) 

DT(2)  =  outdoor  coll  fin  tip  diameter  (In) 

RPCH(2)  =  pitch  between  tubes  of  the  same  depth  In  outdoor  coll  (In) 

DPCH(2)  =  depth  pitch  for  outdoor  coll  tubes  (In) 

WITDTH(2)  =  outdoor  coll  width  (equal  tube  length)  (In) 

Line  34;     FPCH(2),  FTK(2),  FMK(2),  TMK(2),  AMAS(2) 

FPCH(2)  =  outdoor  coll  fin  pitch  (In) 

FTK(2)  =  outdoor  coll  fin  thickness  (In) 

FMK(2)  =  outdoor  coil  fin  material  thermal  conductivity 
(Btu/(ft  •  h  •  F)) 

TMK(2)  =  outdoor  coll  tube  material  thermal  conductivity 
(Btu/(ft  •  h  •  F)) 

AMAS(2)  =  air  mass  flow  rate  through  outdoor  coil  (Ib/h) 

Line  35;     C0NST(2),  CP0W(2),  ANGLE(2) 

C0NST(2)  =  constant  for  air  side  heat  transfer  correlation  for 
outdoor  coll  equal  to  0.134  (-) 

CP0W(2)  =  constant  for  air  side  heat  transfer  correlation  for  outdoor 
coll  equal  to  0.681  (-) 

ANGLE(2)  =  angle  between  outdoor  coll  face  and  air  streamlines  (rad) 
Line  36;  EIDFAN 

EIDFAN  =  outdoor  fan  energy  input  rate  (kW) 
LINE  37;  NREPTO 

NREPTO  =  number  of  repeating  sections  in  outdoor  coil  (-) 
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Line  38;     NTUBE(2,I),  I  =  1,5 

NTUBE(2,1)  =  number  of  tubes  in  first  row  in  each  section  of  outdoor 
coil  (-) 

NTUBE(2,2)  =  number  of  tubes  in  second  row  in  each  section  of  outdoor 
coil  (-) 

NTUBE(2j3)  =  number  of  tubes  in  third  row  in  each  section  of  outdoor 
coil  (-) 

NTUBE(2,4)  =  number  of  tubes  in  fourth  row  in  each  section  of  outdoor 
coil  (-) 

NTUBE(2,5)  =  number  of  tubes  in  fifth  row  in  each  section  of  outdoor 
coil  (-) 

Line  39:     IFR0M(2,I),  I  =  1,10 

IFR0M(2,1)  =  number  of  tube  of  outdoor  coil  from  which  tube  1 
receives  refrigerant  when  outdoor  coil  works  as 
evaporator  (-) 

IFR0M(2,2)  =  number  of  tube  of  outdoor  coil  from  which  tube  2 
receives  refrigerant  when  outdoor  coil  works  as 
evaporator  (-) 

IFR0M(2,3)  = 


IFR0M(2,9)  = 

IFROM(2,10)  =  number  of  tube  of  outdoor  coll  from  which  tube  10 
receives  refrigerant  when  outdoor  coll  works  as 
evaporator  (-) 

Line  40;     IFR0M(2,I),  I  =  11,20 

IFR0M(2,I)  =  number  of  tube  of  outdoor  coil  from  which  tube  I 
receives  refrigerant  when  outdoor  coil  works  as 
evaporator  (-) 


Line 

41; 

IFR0M(2,I) , 

I  = 

21,30 

Line 

42; 

IFR0M(2,I) , 

I  = 

31,40 

Line 

43; 

IFR0M(2,I) , 

I  = 

41,50 

Line 

44; 

IFR0M(2,I) , 

I  = 

51,60 
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Line 

45: 

IFR0M(2,I), 

I  = 

61,70 

Line 

46: 

IFR0M(2,I), 

I  = 

71»80 

Line 

47: 

IFR0M(2,I), 

I  = 

81,90 

Line 

48: 

IFR0M(2,I), 

I  = 

91,100 

Line 

49: 

IFR0M(2,I), 

I  = 

101,110 

Line 

50: 

IFR0M(2,I), 

I  = 

111,120 

Line 

51: 

IFR0M(2,I), 

I  = 

121,130 

IFR0M(2,I),  =  number  of  tube  of  outdoor  coil  from  which  tube  I 
receives  refrigerant  when  outdoor  coil  works  as 
evaporator  (-) 

Line  52:     CAPIDl,  CAPLl,  NCPLl,  CAPID2,  CAPLl ,  NCPL2 

CAPIDl  =  inner  diameter  of  cooling  operation  expansion  device  (in) 

CAPLl  =  length  of  cooling  operation  expansion  device  (in) 

NCPLl  =  number  of  cooling  operation  expansion  devices  (-) 

CAPID2  =  inner  diameter  of  heating  operation  expansion  device  (in) 

CAPL2  =  length  of  heating  operation  expansion  device  (in) 

NCPL2  =  number  of  heating  operation  expansion  devices  (-) 

Line  53:     YL,  YD,  YKl,  YDl,  YK2 ,  YD2 

YL  =  length  of  compressor-outdoor  coil  tubing  (in) 

YD  =  inner  diameter  of  compressor-outdoor  coil  tubing  (in) 

YKl  =  thermal  conductivity  of  compressor-outdoor  coil  tubing  material 
(Btu/(ft  •  h  •  F)) 

YDl  =  outer  diameter  of  compressor-outdoor  coil  tubing  (in) 

YK2  =  thermal  conductivity  of  compressor-outdoor  coil  tubing 
insulation  (Btu/(ft  •  h  •  F)) 

YD2  =  outer  diameter  of  compressor-outdoor  coil  tubing  insulation 
(Btu/(ft  •  h  •  F)) 
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Line  54: 


RL,  RD,  RKl,  RDl,  RK2 ,  RD2 


RL  =  length  of  compressor-indoor  coil  tubing  (in) 

RD  =  inner  diameter  of  compressor-indoor  coil  tubing  (in) 

RKl  =  thermal  conductivity  of  compressor-indoor  coil  tubing 
material  (Btu/(ft  •  h  •  F)) 

RDl  =  outdoor  diameter  of  compressor-indoor  coil  tubing  (in) 

RK2  =  thermal  conductivity  of  compressor-indoor  coil  tubing  insulation 
(Btu/(ft  •  h  •  F)) 

RD2  =  outer  diameter  of  compressor-outdoor  coil  tubing  insulation 
(Btu/(ft  •  h  •  F)) 


Line  55:    RYL,  RYD 


RYL  =  liquid  line  length  (in) 
RYD  =  liquid  line  diameter  (in) 


Note  on  coil  input  data: 


As  a  coil  depth  it  is  understood  direction  perpendicular  to  a  coil  surface 
facing  incoming  air.    Line  13  and  33:     for  fin  tip  diameter  DT(I)  refer  to 
figure  21.    Lines  19-31  and  39-51;  tubes  are  numbered  consecutively  from  the 
first  tube  in  the  first  row  (facing  incoming  air)  to  the  last  tube  in  the  last 
row.     Enter  0  (zero),  if  considered  tube  receives  refrigerant  from  coil  inlet 
port;  enter  999,  if  the  tube  is  nonexistant. 
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Table  15.     Example  of  Heat  Pump  Data  in  the  Input  Format  to  Program 
HP  SIM 


0$0$0S*SYST42( 1 ) 


1  SYSTEM  2  **♦* 

2  5. .5.228,0.744,0. 1743 

3  0.1 996 .0. 3355 , 0 .5446 , 0 . 6540 , 0 . 6958 

4  0. 7233, 0. 74 14, 0.7519, 0.7 595. 0. 762 3,0.  7633 

5  0.05,0. 1 ,0.2.0.3,0.4 

6  0.5,0.6,0.7.0.8.0.9, 1 .0 

7  3578.  , 23  .  53 , -831 . 4 .3542 .  . -6140.  . 3234. 
8'  4.7  73E-7. 1 . 353E -  5 . 2 . 7 55E -  4 . 9 . 436E  -  5 

?  37 . 16 , 105.7 .0 . 231 2 . 0 . 221 3 . 3 .  1  70 

10  0.905. 1 .637E-05 

11  10. ,4. . 100. .0. 

12  3.29 

13  . 3 1 5, . 375 . 1 .055. 1 . ,0. 875 , 1 8 . 

14  .0769 , .0055 , 1 18 ., 118 .. 5800 . 

15  . 1 34. . 681 . .642 

16  .640 

17  1 

10  29,30,29.0,0 

19  2, 32. 4,0, 34. 5, 6, 7. 0,9 

20  10,11.12.13,16,17,18,19,20,21 

21  22,23,24,25.55.0.26.27,58.31 

22  1,3,3,33,34,35,36.37,38.39 

23  40,41  .42.  14.15,47,48,49,50,51 

24  52 , 53 . 54 , 25 . 56 . 57 , 27 . 28 . 29 . 30 

25  60,61,62.63,64,65,66,67,68.69 

26  70. 71 . 43 , 44 ,45 , 46 . 78 . 79 , 80 . 81 

27  82 , 83 , 84 , 85 , 86 , 87 , S8 . 59 . 999 , 999 

28  999,999,999,999,999,999. 999,999,999,999 

29  999 , 999 , 999 , 999 , 999 , 999 . 999 , 999 , B99 , 999 

30  999 , 999, 999 .999 . 999 , 999 . 99S , 999, 999 , 999 

31  999 , 999 , 999 , 999 , 999 . 999 , 999 , 999 , 999  ,  999 

32  2,36 

33  0.315.0.375, 1 . 128. 1 ., 1 . .64. 

34  D. 062 5, 0.0055. 1 18. . 1 18. . 16600. 

35  .134, .681 , .671 

36  .510 

37  6 

38  6,6,0.0,0 

39  2,3.4,5.6,0,1  .7,0,9 

40  10 . 1 1 , 999 , 999 , 999 . 999 , 999 , 999 , 999 . 999 

41  999 . 999, 999 .999 . 999 . 999 , 999 . 999 , 9  99 . 999 

42  999 , 999, 999 ,999 , 999 , 999 . 999 . 999 , 999 , 999 

43  999 ,999 , 999 ,999 ,999 . 999 . 999 . 999 , 999 , 999 

44  999 .999, 999. 999 .999 . 999 . 999 . 999 , 999 . 999 

45  999 , 999 , 999 , 999 , 999 , 999 . 999 , 999 , 999 , 999 

46  999,999,999.999,999.999,999,999,999,999 

47  999.999.999.999.999.999.999,999.999.599 

48  999 , 999 , 999 . 999 , 999 , 999 . 999 , 999 , 999 . 999 

49  999 ,999, 999 , 999 , 999 , 999 . 999 , 999 . 999 . 999 

50  999 . 999. 999 .999 . 999 , 999 , 999 , 999 , 999 , 999 

51  999 , 999 , 999 , 999 , 999 , 999 , 999 . 999 , 999 , 999 

52  0. 078 ,0.5,1 ,0. 073 .0. 5 . 1 

53  50. .0.68,223. ,0.750,0.05. 1 .75 

54  100. ,0.68.223. ,0.750,0  05, 1 .75 

55  180. ,0.311 


END  PRT 
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COP  -  coefficient  of  performance  (-) 

ELUSE  -  total  heat  pump  energy  consumption  rate  (kW) 

QLOAD  -  net  capacity  of  the  heat  pump  (Btu/h) 

RMASS  -  refrigerant  mass  flow  rate  (Ib/h) 

TMASS  -  refrigerant  charge  (lb) 

TG3  -  refrigerant  saturation  temperature  at  compressor  can  inlet 

(F) 

TG6  -  refrigerant  saturation  temperature  at  compressor  delivery 

valve  (F) 

TSUP3  -  refrigerant  superheat  at  compressor  can  inlet  (F) 


The  last  three  output  data  TG3 ,  TG6 ,  and  TSUP3 ,  are  final  values  iterated  from 
TG3,  TG6,  and  TSUP3  supplied  as  guessed  input  values. 

1.3     EXAMPLE  OF  RUN  OF  HEAT  PUMP  SIMULATION  PROGRAM  HPSIM 

The  following  is  a  computer  printout  for  a  HPSIM  run  corresponding  to  laboratory 
test  No.  2.2  of  this  study  (refer  to  table  4). 

Input  in  this  run  was  as  follows : 

-  refrigerant  data  -  constants  for  R-22  as  given  in  table  1.3 

-  heat  pump  data  -  as  given  in  table  1.5. 

-  run  controlling  data,  operating  conditions,  guessed  refrigerant 
parameters  as  shown  on  printout  (lines  1-150,  606) 

The  solution  was  iterated  in  4  energy  balance  loops  (lines  151-252,  253-355, 
356-458,  459-561).     Final  simulation  results  can  be  found  at  lines  570  through 
601.     Computation  of  this  example  required  112  seconds  on  Univac  1108  computer. 
In  the  given  example  there  was  no  need  to  iterate  refrigerant  superheat  at  com- 
pressor inlet.     For  cases,  where  such  iteration  is  required,  iterations 
presented  here  would  be  repeated  3  or  4  times. 
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